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CHAPTER [I 


INTRODUCTION TO DISCRETE ORDINATES METHODS 


IN TRANSPORT THEORY 


Transport Theory 

The behavior of a nuclear reactor and the effective- 
ness of its shield are governed by the distribution in 
space, velocity, and time of the neutrons and photons in the 
system. The transport equation is a statement of conserva- 
tion for these neutral particles, the solution of which 
determines their distribution in phase space. Derivation 
of the transport equation is presented in most transport 


theory texts.l^? 


4,5 but 


A few special systems have exact solutions, 
for more practical systems, numerical solutions of an 
approximate transport equation are sought. A few of the 
more useful transport approximations are spherical harmonics? 


8 and moments methods.” 


Monte Carlo,’ discrete ordinates, 
The spherical harmonics method treats the anisotropy 
of the particle distribution and cross sections to various 
degrees of approximation. Its application to plane, spher- 
ical, and cylindrical geometries has been useful but for 


more complex geometries, the spherical harmonics method is 


SO complicated that other methods are used. 
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The moments method, used in shielding studies, 
calculates particle transport rigorously but the method is 
limited to infinite homogeneous media. 

The discrete ordinates and Monte Carlo methods are 
the most nearly rigorous for multidimensional problems but 


10 ana is currently 


the former is hampered by ray-effects 
limited to two-dimensional problems. The Monte Carlo method 
is the only useful technique for rigorous solution of 


three-dimensional problems. 


Development of Discrete Ordinate Method 
The use of discrete ordinates was first suggested by 


Wick in 194311 and developed for radiation transport in 


12,13 


stellar atmospheres by Chandrasekhar. Carlson intro- 


duced the angular "segmented" Sy approximation for radiation 
transport calculations at Los Alamos in late 1952 and early 
1953.14 A complete history of Sy development is given by 


15 


Carlson and a current bibliography is available.1® 


x17 gives a detailed historical review of reactor 


Chernic 
Physics calculation methods from the Manhattan Project of 
the 1940's to the present and describes the place of Sy 
computer codes in reactor criticality and shielding appli- 


cations. 


Early Sy methods have been refined to the present 


8 which 


“diamond difference" scheme, described by Carlson, 
is incorporated in most Sy codes. New approximations, 


Such as ٣ ۹5 have been derived for special problems. 





a 


Phe L967 eer attempted to derive an Sy differ- 
ence scheme, called VSy, from a variational principle but 
produced no improvement over existing methods. More 
recently, effort has been expended to derive discrete 
ordinate approximations from space-angle synthesis tech- 
niques. 2° Kaplan“! derives an S, approximation in X,Y 
geometry in this manner by using step functions in direction 
as trial functions to eliminate ray effects. Natelson?? 
derives discrete ordinate approximations from a first-order 
variational principle in X,Y geometry. For some unknown 


reason, this approach has not led to improved difference 


schemes over the standard diamond difference Sy method. 


Applications of Sy Method 

Transport theory, as opposed to one of the simpler 
approximation techniques such as diffusion theory, must be 
used when either very precise solutions are required or 
flux variations are large. As a rule of thumb, transport 
theory is used when the logarithmic gradient of the particle 
flow density is of the order of or larger than an inverse 
mean free path over significant portions of the system. 
Early application of SN by Lee23 to the one-speed problem 
revealed that the Ss approximation gave critical radii 
within 0.3% of the values obtained by the exact solution of 
the neutron transport equation. Problems exist, however, 
for which the one-speed approximation is not sufficient. 


24 


Mills applying Sy to small fast critical assemblies such 


as Godiva, Jezebel, and Topsy, revealed that 24-group, Sg 
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calculations with a P, approximation to the scattering 
anisotropy should be sufficient to provide an accurate 
treatment of neutron transport in fast metal assemblies. 

More recently, the demand for tighter safety 
standards has awakened renewed interest in application of 
transport theory Sy approximations. Protsik?? compared 
results of multigroup diffusion theory and Sy calculations 
and found that significant error in loss-of-coolant reac- 
tivity calculations occurs using diffusion theory. He 
found that diffusion theory over-predicts the neutron leak- 
age and under-predicts the multiplication. 

The Sy approximation has also been used as a standard, 
for problems without an exact solution, by which the merits 
of other techniques are judged. 2028 Sy solutions of the 


adjoint equation have been used in importance sampling tech- 


niques for Monte Carlo methods. 2? 


Multigroup Sy Equations 

The steady-state multigroup transport equations can 
be written? 

Lp (r , Q) = Sy (x, Q) + q (r ,Q) (1) 
where the vectors Y and q contain the unknown distributions 


and sources in each energy group, respectively. The opera- 


tors L and S take the following forms 


(Lv(z, $), = 8٠٢۷و (8,ع)‎ + 0 (x) Vg (x8) (2) 


A لا‎ 


Y Jd £(r:8',9' 一 MO E) (3) 
Aj 
gi & 


(Sp (2,2), 


where the subscript g denotes the ath energy group. 
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The Sy approximation discretizes the above equations 


at N discrete angular ordinates, e to produce the set 


(ህ(፲፡8. ) = 59(፪፡ህ مع‎ (4) 
in which 
(Ly(x,2,)), = Q Vis (c, 8) + T A (5) 
and 
^ N ^ ^ 
(Sple, 29), = E DOC pg > RG) gr lor) (6) 
9 m=) 


where w, are weights associated with the angular quadrature 


m 
set {2 and the sum over m approximates the integration 
over Ds in (3). Equations (4), (5), and (6) describe the 
SN approximation to the transport equation. 
Finite-difference techniques using the diamond dif- 


8 have been applied which approximate 


ference approximation 
equation (4) to second order in all independent variables 
in Cartesian geometries and to order Ar*/r? in one- 


31 Applying any method of 


dimensional spherical geometry. 
spatial discretization, the Sy equations can be written in 


the matrix form 


Ly = Sy + q (7) 
where the vectors y and q now contain the unknown flux and 
inhomogeneous source, respectively, in each group and at all 
mesh points, space and angle. The matrices L and SŠ are 


approximations to the operators L and Š respectively. 





Cater Iteration 


The scattering matrix S is split as follows 


uL S S (8) 


u 


where 54,5 and S, represent down scattering, self or 


A ? 
in-group scattering, and up-scattering respectively. The 


iterative technique 


[1-54-5,)9" = S y? * q (9) 


represents the "outer" iteration where y is the iteration 
index. This matrix equation is never explicitly solved in 
puuuspostocodes. The matrix {L-S,-S,} is dense and of such 
large dimension that it is not easily invertible on present 
day computers. Consequently another splitting is made so that 


an inner iteration can be conducted within each energy group. 


Inner Iteration 
The self or in-group scattering matrix is separated 
from the left hand side of (9) to operate on the pth inner 


iterate group flux, 


Jti k+l DESEE 7 
Lava = 559"9 Í 53 (10) 
where the source to the gen group is 
f í*1 j 
نه‎ * (Sav JE a p * qq - (11) 


In this inner iteration scheme, the first term on the right- 


بر 
hand-side of (11) can be evaluated even though y? is not‏ 





ው ገ 


known for groups g and below, since elements of Sy corre- 
sponding to these groups are zero. The matrices La and 
559 in (10) are restrictions of L and S, to the ath group. 
The vector vg contains components of ys?! corresponding to 
the gth group. The matrix La is triangular and thus easily 
invertible for typical difference schemes, so هو‎ can 
be obtained rapidly from ان‎ 
A converged inner iterate solution in each energy 
group completes one outer iteration. Consequently the 
outer iterate equation (9) is never solved but its solution 
is approximated by completing the inner iteration for each 
group a number of times, each with a source, q, calculated 
from previous outer iteration fluxes. For media which do 
not reproduce neutrons, only one outer iteration is required 
but for fissioning media, many outer iterations are required. 
The problems treated here employ the inner iteration 
of equation (10) for plane geometry. The specific equations 


Will be developed later. 


Convergence of Sy Approximation 

Early versions of the discrete ordinates equations 
were used and results compared to those of other methods 
without regard to their formal convergence properties. 
Numerical experiments with problems representing physical 
systems showed that the early (circa late 1950's) time- 


dependent schemes were stable . 32 


Not until 1960 did proof 
of convergence of the steady-state Sy approximation appear 


when Keller22 proved mean-square convergence of the Sy 
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approximation to the analytic solution in plane geometry, 
one-speed, with isotropic scattering for values of 
0 < C < 1 where C is the mean number of secondary neutrons 
per collision. At the same time, Wendrof£?* proved point- 
wise convergence for the same problem for C < l anda 
weighted mean-square convergence for C = l. Keller?° 
immediately thereafter proved pointwise convergence for the 
problem for values of C greater than unity. 

Not until 1968 did convergence proofs appear for 


36-38 proved pointwise 


multidimensional geometries. Madsen 
convergence of the Sy approximation to the time-independent 
one-speed angular segmented transport solution in x,y 
geometry with vacuum and periodic boundary conditions. In 
addition he proved that the spatially differenced approxi- 
mations using the central difference, first-order difference, 
and diamond difference schemes converge pointwise to the 
angular Sy approximation solution in x,y geometry. Madsen”? 
has shown that the Sy approximation in x,y,z geometry con- 


verges pointwise to the exact solution under certain condi- 


tions. 


Acceleration of Convergence 

Current work is being directed toward improved tech- 
niques to accelerate convergence of the Sy algorithm. 
Clifforat0 gives a lucid account of 5۷ running times for 
shielding problems using various generation digital com- 
puters. With present generation IBM/360 series machines, 


very large two-dimensional problems require 10 to 20 hours 
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running time. For deep penetration problems in media with 


: 30 
o°/o® near unity, Reed 


states that non-accelerated Sy 
methods converge so slowly that they are practically use- 
less. 

Several methods have been developed to accelerate 
inner iteration convergence of the Sy algorithm. Carlson 
and Be1141 proposed a scale factor technique which is used 
in present standard codes like ANISN. It uses a system- 
wide neutron conservation principle and iterates until a 


42 found that 


"false" source is arbitrarily small.? Engle 
the scale factor technique accelerates convergence rapidly 
for absorbing media and near source regions for scattering 
media but is slow to converge at points many mean free 

paths away from the source in predominantly scattering media. 
For the latter problems, the scale factor rapidly approaches 
unity (even when the group flux solution is far from con- 
verged) and thereafter does little to accelerate conver- 
gence. Consequently, he proposed a separate scale factor 
for each space interval which exhibited some success in a 


43 proposed an 


one-dimensional 26-group problem. Clancy 
outer iteration scaling technique which, when used with 
inner iteration scaling, accelerates convergence for prob- 
lems characterized by the presence of significant upscat- 
tering. 


The acceleration technique often used in one-dimen- 


sional problems is that of Chebychev, adapted to transport 
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approximations by Hageman. 44 It is not as effective in 
two-dimensional geometries. 


45 and 


The Synthetic Method, first proposed by Kopp 
applied by Gelbard®® has been quite useful in accelerating 
convergence of the Sy algorithm in one and two-dimensional 
problems. This method uses the solution of a diffusion 
theory approximation to accelerate convergence of the 
Iterative Sy algorithm. 

Another effective method of accelerating convergence 
is called coarse mesh rebalancing. An early version was 
suggested by Wachspress^/ for diffusion theory codes. The 
method has been applied to a transport code, TWOTRAN. T8 A 
variational rebalancing scheme was devised by Nakamura ^? 
and applied to the one-dimensional Sy algorithm. 

The most effective acceleration methods proposed to 
date are the synthetic and coarse mesh rebalance methods. 
When successful, they lead to tremendous reductions in com- 
puting time. Unfortunately, Reed 3° discloses, there are 
problems for which the use of these acceleration techniques 
lead to an unstable algorithm. This failure occurs on 
those problems where it is most necessary to accelerate 
convergence, namely problems with an optically thick region 
with scattering ratio near unity. He displays model prob- 
lems for which the unaccelerated So algorithm requires more 
than 2600 iterations to converge and the above two methods 


fail to converge. He also generates parameters which force 


both methods to converge. 





Ze 


MEESbiUeunefor which acceleration techniques‏ ت29۰ 
are required is neutron propagation through large thick-‏ 
nesses of iron which has a huge scattering resonance at‏ 


about 25 kev. Devillers?? 


compares Sy codes running times 
(NIOBE, 20 minutes; ANISN, 8 minutes) to a Monte Carlo code 
(POKER, 13 minutes) for a one-dimensional multigroup prob- 
lem. Deep penetration problems in iron have also been 
peeved by Su techniques to study the "window" effect for 
such materials.>l 

Current research is directed toward developing more 
accurate transport results. Increased emphasis on reactor 
safety is one factor responsible for this trend. But it is 
not sufficient to strive for accuracy alone; efficiency or 
cost must also be considered. Consequently, acceleration 
techniques play a prominent role in present computer codes. 
Any technique which contributes to improved efficiency will 
be incorporated to reduce the expensive computer computa- 


tional costs. 


Proposed Investigation 

This investigation formulates the angular flux form 
of the Sy equations in plane geometry within a mathematical 
frame which allows an insight into the mechanics of the 
convergence process. This mathematical formulation leads 
to sufficiency conditions under which iterative convergence 
to an exact discretized solution is assured. These condi- 
tions include a domain in which some elements of the 


iteration matrix are allowed to be negative. 
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The infinity-norm of the iteration matrix is calcu- 
lated for a problem with a homogeneous isotropically scatter- 
ing medium. This allows the calculation of an improved con- 


vergence criterion. Reed”? states that the convergence 








criterion 
وت‎ 
max " 
2 
where gr is the ¿th component of the kth scalar flux 


iterate is not, strictly speaking, a measure of the error in 
the final iterate, since slowly converging methods may meet 
this test without being close to the exact solution. It 

is, however, the convergence criterion actually used in 

most transport codes. This investigation proposes an 
improved convergence criterion which, along with the 
demonstrated convergence properties of the S, algorithm, 
guarantees that the fractional iterative error, defined in 
Chapter II, is arbitrarily small. This is exhibited in 
Enapter II. 

A spatial transform method which accelerates inner 
iteration convergence for certain S, problems is developed 
in Chapter III. The method is applicable to problems for 
optically thick media in which the scattering ratio is near 
unity. The transform renders the angular discretized Sy 
equations invariant and places an acceleration parameter in 
each non-zero matrix element of the iteration matrix such 
that its norm is reduced. This results in accelerated con- 


vergence of the Sy algorithm. 
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The reciprocal of the spatial transform is applied 
fo an Sy problem to reduce the discretization error. This 
variation of the transform method is applied to a predom- 
inantly absorbing medium problem in Chapter IV. 

An upper bound is found for the round-off error in 
the Sy iterative algorithm describing an inner iteration of 
a multigroup problem. The expression for this bound pro- 
vides insight into the reasons why iterative techniques 
often exhibit smaller roundoff error than competitive matrix 
inversion methods and how these errors propagate through the 
iterative process. This exposition is presented in Chapter 


V. 


mE 





CHAPTER II 


DISCRETE ORDINATES AS A METHOD OF 


SUCCESSIVE APPROXIMATIONS 


To analyze the convergence properties of the Sy 
algorithm, it is necessary to provide some theoretical back- 
ground and to show that the algorithm possesses, under 


certain restrictions, sufficient conditions for convergence. 


Theoretical Background 


The problem to be solved is 


Ly*= q (1) 


where L is a finite dimensional square matrix, y* is the 
solution vector, and q is the inhomogeneous vector. All 
vector and matrix elements are real. For the applications 
considered here, L is non-singular so equation (1) has a 
unique solution, y*, for each given q. 


Now assume that L can be split as follows, 
L = D-E-S . (2) 


The choice of the splitting is arbitrary but the implication 


| but (D-E) is much 


is that (D-E) 7! is an approximation to ا‎ 
Simpler to invert than lL. For instance (D-E) may be chosen 


such that it is sparse compared to the dense L. 


2ے 





E 
After splitting L, equation (1) becomes 
(D-E) [1- (D-E) S]y* = q 
or 
y* = [TI-(D-E) -TS]-1(D-EB) -4 
where 
L^! = [፲-(5-5)” !5]1”!(5-5)”፤'. 


Use will now be made of a convergence theorem from 


53 


Isaacson, which is repeated here. The notation has been 


modified for convenience to conform to that of the Sy 
algorithm used later. 


Theorem 1 


00 


The geometric series ፲[5-51”'5 
m= o 


m ! 
] converges if and 


only if 
[11(D-E) ”Is < 1 (3) 


If (D-E)^!s is convergent, then [7=(D=E) 531 is non-singular 


and 


1۹ ር E FD E) S” 


=5 mm8 


人 
0 
If the splitting of L is such that expression (3) is 


Satisfied, the theorem can be invoked to produce 


Ux! = 15-51 !51"(ሀ-5)”፤ 
m=0 


mE 





عم ےت 


| 


which is the Neumann series expansion of L ` and 
00 
AEO El SI D-E) Tg (4) 
m= o 
A consequence of the series expression for LT) is that 


-1 
٣ lee) dI... (5) 


1-| | sl 
Observe that ||L^!|| is bounded provided that | | (D-E)7'S| | 
« ] and that the bound becomes tighter as || ወ-5)”፤511 


approaches zero. 


Iterative Process 
On the basis of Theorem l, Ra11>* constructs an 
iterative process, called the method of successive approxi- 


mations, which has 


ll 8 (p-E)~!sy(4) + (D-E) lq (6) 


where 4 is an iteration index. This process gives a 
sequence (y €) y Of successive approximations which converge 
ES the exact Solution yp* starting from any initial guess 
plo). Equation (6) implies that 


e 


: Á | 
yet] = ERE DDE lo + [ (D-E) Is glo). (7) 
m=0 


Observe that Lim y (+1) = y* since the first term on the 
L> مه‎ 5m 


right hand side of (7) approaches wy* and the second term 


vanishes due to condition (3). 





uj - 


Error Estimation 
Using equations (4) and (7), it follows that the 


absolute error, defined by 
اك‎ (8) 


is bounded by 


| "ከ1 11 1 ٢ 
IARE = + |] (0-5)7?g]| (9) 


finer ie S| | 
for "1 ES 


This error expression is of little practical value except to 
illustrate the idea of efficiency. It is desirable to 
obtain approximations of a given accuracy with the fewest 
number of iterations. Obviously, more iterations are 
required to achieve a given accuracy as moa 


approaches unity. Indeed the error is unbounded for 


|] -B)7!s]] = ]. But this condition violates expression 
(3) which guarantees that the series converges to Ll. If 
89-51 '5|| > 1, the convergence of the iterative process 


is not guaranteed. 


Formulation of Sy Algorithm 
The Sy equations will be derived for an energy group 
in a multigroup problem applying the following assumptions 
to the transport equation: 
1) steady state 


2) infinite plane geometry with one spatial dimension 
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3) homogeneous medium 


4) isotropic scattering in the Laboratory Coordinate 
System 


5) no reentrant current at the slab boundaries. 


The problem may be described by 


l 
g^ 
ېم و‎ ,  ,٤ؤ,‎ + ql) Q0 


zd 


with boundary conditions 


ylo,u) = O for u >< O (lla) 


u 
° 


v(L,u) for u < O (11b) 


where X = 0 at the left face of the slab and X = L at the 
right face of the slab. Equation (10) is the balance 
equation for the angular flux in a particular energy group. 
The source term q(X,u) represents neutrons scattered from 
other energy groups and born within the energy group, either 
from fission or external sources. 

The angular and spatial discretization of equation 
(10) follows the procedure described by 8611, Alternate 


26 who applies numerical difference 


formulations by Mynatt, 
methods to the analytic balance equation, and Carlson,? who 
immediately writes a discretized form of the balance 
equation with unknown coefficients and evaluates them by 
applying neutron conservation principles, produce the same 


difference equations in plane geometry. Equation (10) is 


discretized in direction cosines by approximating it at each 





OE 


of the N values of uj in the quadrature set selected. This 


gives the N coupled set of equations 


dy .(x) 


4 እ 
12 
u; u 9^ ; (x) - 7 Jupbylx) + ٤17 E 
n=] 


(13) 


for j 1,2,°°°, N where UN SL ee 


; J 
7*j 
for Ebo ec MN Z. 


A spatial mesh is selected with R equi-spaced incre- 
ments providing R+] space points, including the boundaries, 


at which y ; Ux) is to be approximated. Each increment, A, is 


defined by 
A = X - X : (14a) 


The x,u mesh is described in Figure l. 
The quantity ሃ/፤።) in equation (12) is approximated 


at each spatial midpoint of the mesh 


X + X 
Eo ٦ (14b) 
by 
. + . 
Üh+1/2, j 5 Pel, j E (14c) 


and the derivative term is approximated by the central dif- 


ference 


D Pert j O Yki (15) 
X 站 


Xk+1/2 
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Direction set properties: 


1) Hy = "Uç, j s 1,2,***,N/2 ? Hf > 0. 
7*j 
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This approximation of V Ex] by ሃኔ í leads to a type 


53 shows 


of error called discretization error.  Isaacson 
that this approximation is accurate to order A*, O(A2). A 
desirable strategy is to choose A small enough that dis- 
cretization error is reasonably small. This choice is based 
on a guess of how the solution is expected to vary across 
the slab. 


The above substitutions applied to equation (12) 


produce the following set of equations. 
For uj > 0 


4 


aja 


N 
K y ea Yen! * erry, j (162) 


N 
WE - *”ክ 7 LU رل ل ان‎ 95+፲/2,7 (?*) 


k = R,R-1,***,1. In these equations 
juz] + ^o% 
d. E GT (16c) 
1 
and 
AS = Aot 
(16d) 
Í 
and 


tk+1/2,j = 4lXk+1/2, Uj) . (16e) 


Equations (16) could be arranged in the matrix form 





and the solution y* found by inverting L. This procedure 
is not useful, however, because L is usually dense and is 
of such large dimension for practical reactor problems that 
present day computers do not have sufficient storage 
capacity to invert such systems. Instead of the above pro- 
cedure, an iterative technique is used for reactor problems 


which solves (16) in the following manner. 


Bor u. > 0, 


1 
4 N ን 
((+1) [<+] ) g y) [4] 
= = 7 + + 17 
B"..; "í EMM 4 ን "bel, n NL Bay, E 
and for u; < 0, 
(i+1) en 1 po [4] 
0 8 Á+ zs ٨ Á 
0 — و‎ T4 M kel n * Vg, ud * مو نم۹‎ 7 
where & is an iteration index. These equations are the 
basic equations describing the Sy algorithm. 
The iterative solution of this system, "0۲ approx- 


inmates the exact matrix solution y*. The error in this 
approximation is called iteration error. A practical 


Strategy is to iterate until this error is arbitrarily small. 


Calculation Sequence 
The convention of performing the calculation 
Spatially in the direction of neutron flow, as pointed out 


by Bell,” has been observed in equations (16) and (17). 
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Consequently, iterative and round-off errors are attenuated 
rather than amplified as the calculation proceeds. This is 
observed from the S, algorithm in its machine form. That 


ms, for رھ‎ 2፡1) 


; N : 
CELL E [4፥! ] አ زرد‎ (4) 
Y . = Y + ل( س‎ + q (18a) 
pe = ሂ E an 2 meh k+1/2,n k+1/2,j 
and for u; < 0 
E"! ሽዓ ተ Lo ፦ su 1] + q | በይ) 
R+] (d a DUM 0۶1/27 ۷ R+1/2,j 


The quantities Si are less than unity for all y so any 
numerical errors in a particular angular flux value will be 
attenuated in calculating the neighboring angular flux value. 
Amplification of these errors would have occurred had the 


convention not been observed. 


Iterative Matrix Formulation of the S, Algorithm 


Equations (18) with boundary conditions 


Wa qo 0 for u; > 0 (18c) 
Up. 1 ; = 0 for رھ‎ < 0 (18d) 
and an initial guess MM = 0 for all k,j is the Sy machine 


algorithm for the problem described by equations (10) and 
(11). 
Equations (17) are precisely the same Sy algorithm 


expressed in slightly different form. If these equations 
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for an R+1 by N mesh (including boundary points) are written 
so that the solution set is arranged in the following vector 


order 


IM 


ሠጉ 


Ur 
Va 2 
V3.2 


و ېم0 


"8+፲ ዘ/2 
፻ " | Un u/2«1 = 


"፳-1 ዘ/2፦! 


V1,N/2«*1 
YR, N/2+2 
VR-1,N/2+2 


Y1,N/2+2 


| ما" 


which has dimension R:N, then equations (17) describe the 


matrix system 


me = S y V] teg E (20) 
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The matrix (D-E) is block diagonal of the form 


(21) 





where each block is an R by R submatrix. R is the number 
of spatial increments in the mesh. There are N such blocks. 


Each block B; has the lower bi-diagonal form 


ይይ (22) 


ኤኤ 


where / is the index belonging to uj in the direction set 


from which 7 and dj are calculated by equations (16c) and 


(16d). 
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The S matrix also has a regular block structure of 


the form 





where the notation Scp and Sgp distinguishes the form of the 
block submatrix and the subscript denotes the index belonging 
to all elements in a particular block. The S matrix has 
certain block regularities listed below: 

1) The same index appears in each column of blocks. 

2) Each block is an R by R submatrix of S. 

3) There are N blocks in each row and column of 
blocks. 

4) The upper right and lower left quadrants of S are 
composed of blocks, each with form Sy. 


5) The upper left and lower right quadrants of S are 


composed of blocks, each with the form Sgp. 





Jo. 


The block S form is 
BD; 


(24) 
5. ھ2‎ 
where each matrix element has the same value 
g^ 


th 


and W; is the 4 weight in the quadrature set. Observe 


that SBD; is lower bi-diagonal. 
The block Sep, form is 
A 


(26) 
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which will be called upper cross bidiagonal. The matrix 


elements are given by (25). 


Inverting the Matrix (D-E) 

Due to the simple block form of (D-E) and the fact 
that each block is lower bidiagonal, see (21) and 22), its 
inverse is of a simple form. In fact peal is of the 


block form 


(27) 





where each block, B of (21) is inverted individually and 


1 


placed in the same position in DEN that it held in 
(D-E). 


Each block p can easily be calculated from (22) 


using the method described in Wylie.?/ Bo" has the form 





-29- 


] 
di 
31 8 © 
A 
2 | 
e <= ] 
E "ያፎ 
a? a p 
ብ ፈ 
bth pow 55-? ማዜ ek- . g;! 
uh ak-1 R- 2 
E (28) 
To R- 2 R-3 2 Í 
٢ x + EE በቹ! 
: 22 E 
ar d; d; ለከ... 


which is lower triangular with elements given by (16c) and 
(16d). 
Since (D-E)~! and S are now known, equation (20) can 


be operated on by (D-E)^!S to obtain 


yet]! = DE s ህ!*) + (D-E)! ፡ . (29) 


The matrix (D-E)^!s is called the iteration matrix 
and can be determined from (23), (24), (26), (27), and (28). 
(D-E) ^!s can be formed by multiplying (27) onto (23) block 
by block: That is, due to the block diagonal form of 
DE L, DE ۹ is formed by matrix multiplying 57' onto 


each block in the ¿Ah row of blocks in S. That is, 


Fas 
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MDE) Ts - 





where each block is of dimension R by R. There are N by N 
such blocks in (D-E)^!s. The following regularities in the 
block structure of (D-E) "Ts are observed: 

1) Each block in the ¿th row of blocks has 87? as 
S perator where 4 is the index for ۰٠ 

2) Each block in the j^^ column of blocks has either 
RED. or Seo; aS an operator where j is the index on Qj. 

3) The upper-right and lower-left quadrants of 
(D-E)7?s are composed of blocks, each of which includes an 
operator of the form RE 

4) The upper left and lower right quadrants of 
IE) s are composed of blocks, each of which includes an 
operator of the form 289/' 

The matrix form of each block of (p-E) Ís is one of 


two kinds, PiSBD, or BE Cp" 


The form of BCSBD; is 
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where e; and d. are values of equations (16c) and (16d) 


ht in the quadra- 


weig 


evaluated at P and 7 is the ¡e 


ture set. 
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(32) 


BS 


= ad دك‎ 








The fact that there are, at most, two terms in each 
element of (p-z)~'s is attributed to the form of DE and 
EU. each of which have, at most, two non-zero elements in 
each column. These properties give an ordered regularity 


to each block of lp-E! 15 


Equivalence of Sy Algorithm to Matrix Formulation 

As far as the author can determine, the Sy algorithm 
has not previously been described explicitly in the matrix 
form of the preceding section. The discovery that the Su 
algorithm, described by equations (18) is equivalent to 
solving matrix equation (29) at each iterate step should be 
no surprise since equations (18) were derived directly 
from equations (17). At each iterative step the Sy algo- 
rithm provides a recipe for inverting (D-E) and multiplying 
(4 } 


the result into Sy and 4. 
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This discovery, however, allows the Sy algorithm to 
be classified as a method of successive approximation which 
provides a mathematical description of its convergence 
properties and error analysis. That is, equation (29) is 
precisely the same as (6). Consequently, properties 
previously discussed for the method of successive approxima- 
Bons apply to the Sy method if the conditions of Theorem 1 
are met. Before this is attempted, it is necessary to show 
that the norm of the iteration matrix can be calculated for 


the Sy problem previously described. 


Calculation of sl 
The norm imposed upon the system is the infinity 


norm, defined for a vector by 


Ill. = maxlypl . (33) 
R 
Obviously, this norm is conveniently chosen for pointwise 
error analysis. Isaacson?? shows that for a matrix, Á, 
||ላ||ሬ = max lag yl , (34) 
ፌ J 


that is, the maximum row sum of absolute values of the 
matrix elements. From this point, the above definitions 
will hold for the norms simply denoted by || ||. 

The row sum of (D-E)"I!s, described by (30), (31), 
and (32) will be a row sum in a row of blocks which com- 
prise (D-E) Is, Since each row of blocks has the same 


form, it is sufficient to sum a row of blocks and pick the 





Ea 


row from the block sum which gives the maximum row sum of 
absolute values for (D-E]-!S. 
Due to the ordering of the direction set, 
UN = W; (35) 


J 
—+ 
7 J 


I - 1,2,---,N/?7. A property of the direction set is 


that 

/ 

; للا‎ 一 ] ° (36) 
Consequently a block row sum of (D-E) 71'S can be collapsed 
by taking 


I -1 -] 
5: 589, + 5 SBD, SA By RD = 


9 - 1 - 1 - ] -] - -] 
(B; 355, * 8; Seo, )*(8. $gp, * 8; Sep, (+. ..+)8 358p, B coh 
m qi 7 


E 
EUG =] 2 =] - -] = 
x ሸገ تہ‎ SEIT BT a Sl O 05 Sop) 


2 
- 9^ (y gos. LH) B2 Se BT I Sn] 
4 MZ up COED CCD 
2 
: Aine Ge 
4 < B) < “CD (37) 


where B: !$gp 1 B;!Scp = sum of the matrix parts of equations 
(31) and (32). 

It is observed from (37), (31), and (32) that the sum 
of absolute values of the elements of each successive row 


include all elements of the preceding row plus a positive 
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quantity. This is true of all row sums except the last. 


The last row sum is equal to 












































- | - 2 
2 (|: : e; D 1 % + ci ..... 1ሟ ) m 
q ái -R “፳-1 . 
The next to last row sum is equal to 
Re kan 
9^ 5 = + 4 1 ۹ “ፈ + ez +° ° ° + “ፈ 
4 aR-1 d; al 3 _R-2 (39) 
Á 7 3 an : 
The last row sum is the maximum if, 
- 2 _ R-2 R- 1 
ei 3 | > Il or if E 3 ^ E (40) 
ا‎ ak an 3 d; 
ፈ. | 











t 
Even if A,o , and |u; | were adjusted so that €; < 0, the 


quotient lec is less than unity because m d. for all. 
. ፍች for systems with large R, the next to last 
row sum is the maximum. 

Now the block row sum must be picked from all possi- 
bilities. The index « for which expressions (38) and (39) 


are maximum can be determined from equations (16c) and (164). 


The quantity 


£ 
Ja | E ? H; 
d; ] + 
4 
attains its smallest value for a given A and o* when Im 


is the minimum in the quadrature set (u;). Also 
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. The 





] ZA 

Q0 0100070777774 iS a maximum for the same |y; 
latter factor dominates, however, in the norm determination. 
Therefore the norm is calculated by using the minimum 


absolute ordinate value in the set tug). Consequently 



































R-1 2 - E 
Da 3 1 + = +4 SE + 4 el I -4ለ ብ ቀ3 | لح«‎ 
Tid} a, 4:۱۱١۱ fa; d; d; qj 
- ] 
Iito-E} SIII < 
(41) 
R-2 2 R-3 

A e. e. e. e. ; 

0 3 | 4 4 ge 7)۹ .4. otherwise 

q B d; i ፦=! al lal پ3ا''۔.'‎ ) 


























where & is the index for the minimum absolute ordinate 
Value in the set tuj). The latter value is appropriate 
r large K. Observe that equality occurs when e;»o. 

For a given A,0%,0%, luslméne and R the infinity- 
norm of the iteration matrix can be calculated by (41). It 
is obvious that the norm is bounded since = < 1 and = 
1s bounded for any finite dimensional system. 


x 


Observe that, for a fixed A,0*, and BEL the 


coefficients 3 and š are fixed. In this case the norm 
decreases linearly with ር Consequently the norm of the 
iteration matrix is small for problems involving primarily 
absorbing media and large near unity for predominantly 
scattering media. This results in rapid convergence to an 
acceptable iterative error for absorbing media and slow 
convergence to the same level of error for predominantly 


scattering media. This behavior is predicted by equation 


(9). It will be verified experimentally later. 
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Sufficient Conditions for Convergence of Sy 

Equation (41) will be used to provide sufficient 
conditions to guarantee that ٣۶۶۶۰٣٣٦٣ < 1 in order that 
the Sy algorithm can meet the conditions of Theorem l. For 
large R, "۶ ۹ٰ٦ is evaluated by the second expression 


in (41), which can be rewritten 


A R-] 
|| 5-5)”!8|| «ና «(5 | NE (42) 
4d; | M-8 
where 0 > 86 = ا‎ < l پټ‎ (43) 
ፈ 


For large R the terms gR- ! and gR-2 are insignificant and 


could be neglected. But since 


imposing the condition 


or [ 1 ) 7 (44) 
di [=8 

will guarantee that | | [5-5፪)”፤5| | « ] without seriously 
loosening the resulting bounds on A. 


From (16c), (16d), and (43), expression (44) becomes 


d, > g^ + 6,۱ 


or 


x x 


፡|ህ4 | ዘፈክ + 40 > Ao? + 8. - Ag (45) 





E 


This inequality provides the sufficiency conditions for 


convergence of the Sy algorithm. 


= and 


Suppose the problem is defined so that o?,o 


Bel min have been selected and it is desired to impose a 


sufficient condition on the choice of A so that (45) is 





Satisfied. Then the condition || ۳ | < 1 is satis- 


fied. 


t So that €e: < 0. 


Suppose further that Ac 1 


d ከ0 


Then (45) becomes 


1 £ 
A ር ۳۴ 7ለ75 + ለሆ" - 21|፡ቱ,/1 ዘፈክ 


or 


٣‏ ھ2 
. — 
O‏ 


(46) 


sS 212955965 a limit on the negativity of the elements e,. 


x 


If Ao” had been chosen equal to 2|y;| Then 0ق‎ 


min * 
4 

in (42) and | < ] is a sufficient condition for 
Á 

||(p-E) !s|| < ]. But this leads to the inequality 


ለ < |" [min But A = 2{ vil min 
£ 


A 


o 0) O 
| | ፅ t 
so the inequality becomes -gz < — or o < 0”. Conse- 
y ርን 
quently ||(D-E)^'S|| « ! in this case if 
0S9 «1, (47) 


c 
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Finally, suppose Ao و ٹ5ت"‎ is selected. This 
choice guarantees that ej and dy are positive for all y, 
hence (D-E)-!S and (D-E)-*! are positive operators.  Expres- 


sion (45) is then 


1፣ከን።ን ግ Nos ez e mew د‎ dom 


Or 


4 
0 
< ] (48) 
ot 
Consequently, || (pD-E)~/s] | < 1 for all physical problems 


with this choice of A. 


Properties of (D-E) ` 

The matrix (D-E) exists in the construction given 
by (21) and (22). It is a linear operator due to its con- 
struction from the discretized linear equations (17). It 


is a bounded operator since 


۰۰" - 11 


4 


which is bounded for all finite values of o^,o ,lu;lai,, 


and A. 
The matrix Gei. exists if det(D-E) is non-zero. 


From equations (21) and (22) 


Hee (DHE) = 111 رر‎ 


and since 00 > 0 for all j, det(D-E) * QO. 
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A conditions Of Theorem | are met by the Sy 
algorithm if the sufficiency conditions listed in Table I 
are met by the problem parameters. If these conditions 
are met, || 5-፪]”!5|| < 1 and the Sy algorithm converges to 


Ņ*, the exact solution of the discretized set of equations. 


TABLE I 


SUFFICIENCY CONDITIONS FOR CONVERGENCE OF Sy 


Bound on A |Values of o°/ot 
to Meet for Which 
Possible Choice| Value of Sufficiency | Sp Algorithm 
For A e Condition Converges 


Negative 


Positive 





The above table can be summarized by stating that 
the sufficiency condition is met for any A provided that 
cm > ot or Bc whichever is smaller. 

The consequences of the Sy algorithm meeting the 
sufficiency conditions for convergence are that expressions 
(1) through (9) can be used to study its convergence and 


iteration error properties. 
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Positivity 

Positivity of the solution of the Sy equations for 
positive source and non-negative boundary conditions is 
desirable because it is a numerical approximation to a non- 
negative physical quantity. In practice the Sy algorithm 


8 One-dimensional discrete 


has produced negative solutions. 
ordinates codes such as ANISN>® and DTF-IV°? have flux 

fix-up routines incorporated so that negative solution 

values are not allowed to occur. This strategy, although 

of practical merit, is not mathematically rigorous. The 
effect of these arbitrary fix-ups on the iteration error 

is not clear. Lathrop? investigated various two-dimensional 
difference schemes and found that the quest for positivity 

is made at the expense of accuracy and efficiency of con- 


vergence. 861177 


points out that the Sy difference equations 
may lead to negative solution values since they do not 
correspond to a positive operator. 

The preceding section revealed that sufficiency 
conditions for convergence of the Sy algorithm can be met 


which allow the iteration matrix to possess negative matrix 


elements. But convergence to the exact solution cannot be 


guaranteed if A is chosen so that A > 1 Conse- 
quently, there is a limit on the negativity of some of the 
elements of (D-E) Is if the theory is to be invoked to 
guarantee convergence. 

For computational economy reasons it is desirable 


to choose A as large as possible. For the problem 





202. 


previously described, however, A < 20ء‎ imposes an 
absolute upper bound on A, above prone ete sufficient condi- 
tion for convergence is not met. If this bound on A is 
observed, it 1S reasonable to assume that the approximate 
solution لكان‎ will be non-negative if X is large enough 

even though initial iterates may have some negative com- 
Penents. This is because y* is expected to be a non- 
negative vector due to its approximation of ሃ/!።) which is 
non-negative by its physical interpretation and y UE! 
approaches y* as 4 becomes large. This is true even though 
the theory of positive operators, which is used successfully 


47 is not applicable in 


in multi-group diffusion theory, 
this instance. Of course no guarantee exists that y=, the 
discrete approximation to ALS is non-negative. Numerical 
experiments by the author confirm that it is for uniform 
and first collision sources when the sufficiency conditions 
are met. 

If A is chosen so that A < اا‎ however, 
(D-E)-!S and (p-E) 7! are ےت‎ RIS operators 
[4] 


which guarantee a non-negative approximation y for all «. 


Necessary Conditions 

The previously discussed conditions are not necessary 
for the Sy iterative process to terminate based upon 
fractional differences reaching an arbitrarily small value. 
Indeed termination of the Sy algorithm iterative process 


has been observed by this investigator for choices of 





- 43- 


A >> psi for which DEI Sj] >> 1. For some sources, 
e.g. ھ‎ orn source, convergence to a positive solution 
vector occurred. For other sources, e.g. first collision 
source, the solution vector contained negative components. 

It is doubtful whether this process can accurately be called 
convergence since the nature of the resulting solution is 


not clearly defined. 

Even if A is chosen so that A > ووو را‎ in which 
Case all ej are negative, Ip-E)^!s — even and odd 
powers of e¿/d¿. Consequently, the iteration matrix is not 
wholly negative. In this instance the Sy algorithm may or 
may not converge to a non-positive solution depending upon 
the source distribution. For these cases, where some 
elements of the solution are negative, doubt exists as to 
the nature of the solution. 

It is not clear that necessary conditions on A exist 
for convergence of the Sy algorithm. It is clear, however, 
that a solution exists if the sufficiency conditions are 
met and that these conditions allow a limited amount of 


negativity for the iteration matrix. For the numerical 


problems previously cited, these solutions are positive. 
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Effect of A^ on ||(D- E) ^!S|| FIGURE 2.2 


NORM BEHAVIOR WITH A 


Flgure 2.2 reveals the 


effect of A on || (D-E)7?S|| : 


for the following problem 









arameters: ` Region 
—2. where all 
5 ري‎ e; < 0 
Slab width = 10.0 cm “ ፈ 
t 1 Ca tive so- 
= = lutions 
z 2n Region of 
o^ = 0.9 cm Í allowed 






negativity 
|Hi|min = 0.23862 


As R (the number of spatial 


1 : 2 
: A, CA 
increments) decreases, A 
becomes larger and ||(D-E)7?S|| increases. In the region 
where A meets the sufficiency conditions, | | (5-ጀ]”?5|| « ] 


and convergence of Sp to the exact matrix solution is assured. 


A small part of this region allows negative elements in 


(D-E) TS. 
In the region where "4 "٥٦ > 1, termination of 
the iterative process may occur. Solution vectors may be 


positive or partially negative depending upon the source 
distribution. This region should be avoided for practical 
problems unless other evidence is available that the solu- 


tion has validity. 


Effect of Quadrature Set on || |[5-51”!5|| 
Expression (9) reveals that the iterative error is 
less tightly bounded as MESI increases towards 


unity. This means that more iterations are required to 
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reach a specified accuracy. Equation (41) shows that 

|| (571) 7! s] | increases for fixed ሰ as the ordinate 
|!i|mín decreases. This is because = grows larger and 

le, I/a; diminishes as A increases but 1/d; predominates in 
(41) so ||(D-E)"!S|| increases. All quadrature sets have 
the property that | i | min approaches zero as the order N of 
Sy is increased. Consequently, the effect of choosing a 
higher order Sy is to increase ||(D-E]^!S||, hence requiring 
a larger number of iterations to convergence. This loss of 
efficiency might be observed in going from S5 to S, where 
the spread in |Mi|min values is significant but is of no 
practical importance for higher order Sy where |y; |min 


values are quite close. 


Convergence Properties of Sy Algorithm 


The following discussion applies to the Sy algorithm 


when A < 2lilmin, Then 7 and ۰ are positive for all j 
O 
and [D-E)7?S and (D-E) 7Î are non-negative. The sufficient 


condition for convergence is met for all physically 
realizable problems. The initial guess, plo! = 0, is a con- 
venient reference so that comparison of convergence 
efficiency can be made with another method later. 

The Sy algorithm, under these conditions, converges 
in such a way that each vector component approaches its 
exact value in a monotone non-negative manner. This can be 
ascertained from applying equation (7) to the difference 


defined by 





> 4 6፦ 


s tT) : y t1) _ "a | ምን 
Then 
 — ፻[(9-8)”!31"(5-81”! : ^ Lp-E) 7 !s1" (5-5) 7! q 
نی‎ m= o 
= ل‎ (D-E) lq EE 


which is positive for q > 0. Consequently, some positive 

quantity is added to each component of y UC! to produce the 
respective component of 00۸ ۲ 

From the definition of ALS given by (8) and y* given 


by (4), 


| | ር in] 
e yop ll. TLp-g17! S1" (p-g]^!a- Yt(p-s)7!s]" (p-]^!q 
m= 0 m=0 


- Fip-g) ^! s] ^(1*(p-E) ^" s«[(p-E)7!s]2«--- 3 (p-g)7!q. (51) 


which is positive for q » 0. But 


o [ED O gr (52)‏ ۔ (اخغاع اغا 


Hence 
g (447) < eli) (53) 


due to positivity of e UC! and ርፌ 1). 566 ۶19076 2.3 for a 


—— 


convergence diagram. 





ag 


The largest absolute error in each successive iterate 


1S monotone decreasing. From equations (51) and (4), 
lé » [[5-5]”!5]" ሄ* . (54) 

Consequently, 

[4] 


This means that 


Em mT for all k which 
. ፌቶ : : ; (4) 
does not imply that ዬኔ 1s some fixed fraction of e, 
for all R even though (52) ensures that Pa < ا ا‎ for 


all k. That is, it cannot be shown that fit « B d 
for all i,k where ይ is a positive constant less than unity. 


The most that can be inferred from (55) is 


WA eee aee (56) 





This means that the maximum error at each iterative step 
diminishes if || መመዌ1”፤51| < ]. This condition is met for 
convergence. 


Similarly it can be shown from equation (50) that 


0۸+0 1 as na ٠ ህህ (57) 








using the same arguments and with the same consequences 


discussed above for the error. 
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It is not of much practical value to know that the 
error at each iterative step is less than the maximum error 
al the previous step. It would be of more practical value 
to know that the error is bounded by some fraction of the 
difference at each step because the differences can be 


calculated. Equation (51) is 


elé) = [(p-E) روا‎ (I+ (D-E) IS+[(p-E) !s)?2+...)YID-E) lq 


(፲+[9-ጀ]”8+[(5-8]”!5]፤ፃ›፦)[(5-5)?!5]“ 5-5)” 


because 7)۶ ھ۹۹۳٢‎ commutes with each term in the braces. 
So 


o 
due to (50). Consequently, 


ቀ] 
A 2 He Md y] < Jas || 


EE a )58( 
۶ ۰۱۰۱۹ ۶ٔ, Si] 


for all k. This does not provide enough information to 


assure that e 


7 is some fixed fraction of an for all k 
but it is enough to provide a practical link between the 


error and difference which can be exploited to derive a 


useful convergence criterion. 
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Convergence of the Sy FIGURE 2.3 


algorithm is pictorially Sy CONVERGENCE DIAGRAM 





described in Figure 2.3 for ہے سے ےت‎ " 

T 
each vector component. The yan En - 
iterative solution proceeds ۷ 3 
progressively towards the E | ^ J 
exact solution in a monotone : 
positive way governed by 
properties (53), (56), (57), 
and (58). ۷ 


Conventional Pointwise Convergence Criterion 
The pointwise convergence criterion based on 


fractional differences, 


(6*1) (X) 


AS‏ د 


€ 
(i) (59) 
n 


where € is chosen arbitrarily small, is used universally in 
practice for shielding problems. It is especially useful in 
deep penetration problems where y* has elements of very 
small magnitude since it weights these regions more heavily. 
Ihat is, iteration proceeds until the deep fluxes, which 

are the last to meet the convergence criterion, are con- 
verged. Its utility is based on the practical idea that, 

as the differences become smaller with successive iterations, 
the fractional amount added to each solution component is 


smaller than the desired accuracy after some i. Consequently, 
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further iteration will not add significantly to the solu- 


tion. 


It must be ascertained, however, whether this con- 
vergence criterion, if met, implies that the pointwise 
fractional error is sufficiently bounded. After all, it is 
desired to iterate until assured that the pointwise iteration 
fractional error is arbitrarily small, not just until the 
pointwise fractional differences are bounded. 

Suppose iteration proceeds until condition (59) is 


met. Since 


eh 2 elél eliti) 


this implies that 


Ey 
ne) 1S el“! so the above implies 
Mn E 
Or 


ei M ps + [5-፪)-!5 elil (60) 
But by definition 


[(p-E) 7! se 4! j, < 1۱ا جوا ۔(ھ-8)]|‎ 


« ||(ውጩ?!5|| ||ዩ“?|| . 
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Consequently (60) implies that 




















e) « e vil ۱٦د ۴ء‎ le 
for all k. But from the definition 
eU - y*-y 1), the above inequality becomes 
ed! « c Dpg-e 1) * [0-517] ||| ]e €! | 
or 
TTS < epp + E sje] 


for all kR. But if this inequality holds for all k, it 
must hold for j, the value of k which corresponds to the 


maximum absolute value of a Hence 


E 





|፣›5-|| -£17!sll] < e yj” 


or 


JA TA E 


| — (61) 
1 fre - Cl 


where j is the index for which A attains its maximum. 
Observe that (61) does not imply that the fractional error 
is bounded for all components of y*. However, if it is 


arbitrarily assumed that 
* » A 
Ta 7 (62) 


for all k, expression (61) implies that 


SGI 


ET. € 
Vp j+e-||(D-E)-/78|| 
or that 
(i) 
E € (63) 


D p+e- |] (D-5)-18]| 


for all k. But even under the dubious assumption (62), this 


does not imply a satisfactory bound on the fractional error. 


Figure 2.4 displays a plot of the bound versus | DE) Tsj]. 
Observe that, as eR Si] approaches unity, the bound 
FIGURE 2.4 


grows to l. That is, up to 
BOUND ON FRACTIONAL ERROR 
100 percent fractional error 


is implied by invoking (59) as I 
a convergence criterion. This 

criterion is useful and ሯ 
practical for || +29 ۶۷3۸2۴۴ < 5 
0.9 but for larger values it is “E 
a poor choice because the B 
resulting bound on the 

fractional error is too large. s 

° || (0-3 *sl| 1 


An Improved Pointwise Convergence Criterion 

Implicit in the concept of a practical convergence 
criterion is that iteration is terminated only after the 
fractional error is arbitrarily small. Since the exact 
solution is not known, it is necessary to iterate on 
fractional differences, which are known. 11 ٠٦ can 


be calculated for the stated Sy problem so it is possible 
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to establish a convergence criterion based on fractional 
differences which implies that the fractional errors are 
arbitrarily small. 


Suppose iteration proceeds until 


۰۷۷ 2 
En: ie [1-11 -መ)”'51||| (64) 
k 


for all k where € is some arbitrarily small positive number. 


This implies that 


||, («+፤) || € ya 
ME for all k. 
1-|| (D-B}~'s| | 
But repeating (58) 
ተ] 
[4] 1 
o Sn so the above 


8 ۹۲ 


inequality implies that 
(4) 


¿k « 3m VIt) for all h. 


But by definition ህይ) = w* - cs so the above inequality 


implies that 


Or 
elél ce ys 
or finally 
vp (65) 





መ መመ 


for all k. Consequently, upon invoking convergence 
Criterion (64), the convergence properties of the Sy algo- 
rithm assure that the fractional iteration error is 
arbitrarily small. This is precisely the goal for a machine 
convergence criterion. 

For the Sy problem discussed here, equation (41) 
allows the exact computation of || |9-ቺ]”!5|| and € is 
arbitrarily chosen. Of course for problems in which 
DE) |5|| is not obtainable, it is necessary to use (59) 
as the convergence criterion but recognition of its inability 
to guarantee a satisfactory bound on the fractional error 


for problems involving strongly scattering media is impera- 
tive. FIGURE 2.5 


: BEHAVIOR OF IMPROVED CONVERGENCE 
Figure 2.5 reveals that CRITERION 


the convergence constant, 
€ 


ومس 


z 4- 
سا‎ [1-1 | (D-E) 'sı]] becomes : 
E ` = 
vanishingly small as ید‎ 
Gl 
|۱۶٢ ھ٦۱‎ approaches unity. 2 
For predominantly scattering ZH 
ہے۔‎ 
media || [D-E)7?S] is quite 7 
| 
Close to unity. It is neces- 9 
sary to use the improved con- E: 
0 |أ5” 2ع -ه)||‎ 1 


vergence criterion (64) for 

these problems. For weakly scattering media, i.e. o^ /o* E 
0.9, the conventional convergence criterion is acceptable. 
Table II shows that the effect of using the new convergence 
criterion is small for ۔8)]||‎ 111 < 0.9 but overpowering 


for || (p- E) !S|| » 0.9. 
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TABLE II 


IMPROVED CONVERGENCE CRITERION 
FOR g = 1074 


— ې سحې ہج‎ A 
سي سی وو‎ et w t. -[.[ [Í[ ° s< swmsss—=ssrrisik  -.4-. vrs— rV s.s = 


ፐጭ [!-||(5-517፤5| |] 





||15-51”!5|| 





0.5 x 1074 
1079 
10 
1077 


The improved convergence criterion imposes a much 
more severe condition for convergence than does the con- 
ventional convergence criterion, especially for predominantly 
Scattering media. For these problems, convergence is slow 
even using the conventional criterion.  Imposing the 
improved convergence criterion makes a bad situation worse. 
It is necessary, however, if assurance is required that the 


fractional error is arbitrarily small. 


Computer Experiment 

A sample Sy problem, described by equations (18), 
was run on a 360/67 IBM computer with the following param- 
eters: 

Slab width = 10.0 cm 

R = 30 spatial intervals 


ot 1 


ዘ 


1.0 cm 
o4 = varied 


e = 1074 





-56- 


N = 6 ordinates using Gauss quadrature sett 


Uniform source, qj (x! SI GO bor all j.: 
It is worthy of note that any of the more recently derived 


23,62 could have been used but 


mechanical quadrature sets, 
the Gauss set is satisfactory to demonstrate convergence 
properties of the Sy algorithm. This set will be used 
throughout. 

A copy of the computer program and samples of the 
computer output are listed in Appendix A. Table III lists 
the results. The number of iterations to convergence are 
tabulated for various o4/o% ratios using each of three 
different convergence criterion. 

Observe that, regardless of the convergence criterion 
used, problems involving predominantly absorbing media 
required fewer iterations to converge than did those for 
scattering media. The number of iterations to convergence 
is a monotone increasing function of 25 /o* , 

Comparing columns two and three shows the effect of 
converging on the maximum difference instead of the point- 
wise difference. As expected, the maximum difference 
criterion is more severe since it takes more iterations to 
meet it than the pointwise difference criterion for each 


value of ሮ* /ሮ" . This is expected because 


; [<+] ) 
Jo lett) | > E 


A A 
y; ን 


for all h by definition. 
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Comparison of columns two and four shows the effect 
of the improved convergence criterion which guarantees 
that the fractional iteration error is less than €. Its 
effect includes the effect discussed in the preceding para- 
graph and has the additional effect of ከ 2 ۳۲ Both 


effects increase the number of iterations to convergence. 


The effect of ||{p-z)~/S|| becomes more dominant as o4/ot 
increases above 0.9. This is due to the direct linear 
relation 0% has on ۸۸" ۹9۹۹۷۳۳ ۰۳۲ 


A comparison of solutions, one of which met con- 
vergence criterion (59), the other meeting the improved 
convergence criterion (64), was made for the values of 


0% /0* listed in Table IV. 


TABLE IV 


COMPARISON OF u VALUES 
e = 10” 


C 


٢ Agreement of Conventional Convergence Criterion So- 
c^ /g lution with Improved Convergence Criterion Solution 


Agreement i significant figure 


All agree significant figure, many to 3rd 
All agree significant figure, some to 3rd 
All agree significant figure, some to 3rd 


All agree significant figure, few to 3rd 





For this problem ||(D-E)~/s|| = ى/ فى‎ 
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The above table reveals that agreement exists only 
to the second significant figure for oS sot > 0.9. From 
previous discussion it was established that a more accurate 
solution results when a larger number of iterations are 
made. Consequently the improved convergence criterion solu- 
tion must be closer to the exact solution than the con- 
ventional convergence criterion solution because more 
iterations were performed for the former. The above table 
shows that the conventional convergence criterion solution 
agreed to only the second significant figure with the more 
accurate improved convergence criterion solution. This 
illustrates that the expected accuracy of € = 1074 was not 
achieved for large values of 04 /ot using the conventional 
convergence criterion. 

Tables III and IV provide vivid illustrations of the 
necessity for using the improved convergence criterion for 
problems with ۶ S| | > 0.9 if an arbitrarily small 
fractional iterative error is required. 

The overall effect of using the improved convergence 
Criterion is insignificant for problems involving strongly 
absorbing media for which || 5-51]”'5|| is small. The effect 
is marked, however, for predominantly scattering systems 
and results in additional computational work required to 
assure arbitrarily small fractional iterative errors. Con- 
Sequently any scheme which accelerates convergence of the 
Sy algorithm by decreasing ||ወው-5)”፤51|| would be eminently 


useful when invoking the new convergence criterion, especial- 


ly for strong scattering media problems. 





CHAPTER III 


ACCELERATION OF Sy CONVERGENCE BY SPATIAL TRANSFORMS 


General Discussion 
The previous chapter discussed how the Sy algorithm 


converges more rapidly for absorbing than for scattering 


4 x 


media. As G0 decreases, holding o“ constant, the number of 
iterations to convergence decreases regardless of the choice 
of convergence criterion. Based on this physical argument, 
it is expected that a mathematical device which adds 
absorption to the problem will accelerate convergence. From 
a mathematical point of view, an operation on the Sy 
equations which decreases the norm of the iteration matrix 


will improve the efficiency of the algorithm. A new method 


is presented here which incorporates both ideas. 


General Spatial Transform 


A transform of the type 
p(x) = ዕ[2])4 (x) (1) 


has been used to provide an effective importance sampling 
device in the solution of the transport equation by Monte 


Carlo methods. Š 3 


It biases the sampling distribution 
towards more important directions, thereby reducing the 
Variance in Monte Carlo methods. This leads to improved 


Statistical accuracy for a given calculational effort. 
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Transforms of this type have been used by the author 
to accelerate convergence of the Sy algorithm for certain 


problems.  Transforms which have proven successful are: 


ax 
ሃ/!።] ። ቆ/(12]6 › (2) 
W(x) = Ole) anh (Ê) , (3) 
and 
+b = 
ሃ/!።) = pled Inr] (4) 


where Q is an arbitrarily chosen acceleration parameter. W 
and b are parameters chosen to restrict the range of the 
arguments. These transforms are applied to the angular 
discretized transport equations. The most effective trans- 


form of those listed is (2). 


The Half-Range Problem 
The spatial transform accelerates convergence of the 
Sy algorithm for the deep penetration infinite slab problem 
under the following assumptions: 
i) steady state 
ii) homogeneous medium 
iii) angular flux has no azimuthal dependence 


iv) high energy particles are scattered only in forward 
directions 


v) forward scattering is equally distributed in 
direction. 


A consequence of these assumptions is that the trans- 


port equation has the form 
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| 
Lo = nd + all (5) 
0 


valid for 


QUE ME 


The anisotropic forward scattering assumption is valid for 
high energy neutrons and gammas. The assumption that the 
forward scattered particles are equally distributed in 
angle is a simplification of the usual Legendre Polynomial 
approximation and is made strictly for convenience. The 
more accurate description of the forward scattering aniso- 
tropy could be treated by the method. 


The scattering distribution function, under assumption 


V) is normalized by recognizing that {(u'*u) = const and 
1 
[const du = | . (6) 
0 


This simply states that the probability that a scattered 
neutron emerge into some allowed direction is unity. A con- 


Sequence of this normalization is that 
4lu'>u) = 1 . (7) 


Angular approximation is accomplished by discretizing 
equation (5) in the N/2 forward directions of the previously 


described quadrature set. This gives the Sy approximation 


dy .(x 


(x) d 
D aan + کس‎ AN nn + ALS (8) 
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with j = 1,2,***,N/2. Quadrature weights obey 
/ 
? s 1 (9) 


and all H; > 0. The angular quadrature set is the positive 
half of the one used for the full-range problem. 
Spatial discretization is accomplished in the same 


manner described in the previous chapter to give 


^ N/2 
E . 2% 
est, j ej 2 E : 1b1/2,j u 
for j = 1,2,°°:,N/2 and k = 1],2,°°:,R 
where 
٢ IE Ao 
pee m (10b) 
2A 
and 
2M , - Ao 


Es — (10c) 


and R is the number of equi-spaced spatial increments, A. 
The Sy algorithm, as before, solves this system of 


equations iteratively as described by 


(Ger) Ka ©: M (i) 7 
d. EN .-. EM — X) 
MIS j ej, ر‎ 2 o ar, ې‎ 1be1/2,j (11) 
where 4 is an iteration index. This set of equations is 


divided by a to give the computer-solved S, algorithm 





ED 


EN. N/2 : i 
EDU S ا اک فول_1, (1+ء)‎ 
1 - d; | a n) 6 ጨጨ 


Numerical results for this algorithm with the first 


collision source 


5 | ah, (13) 


Ak+1/2,j ت‎ 


and vacuum boundary conditions 


Yi = 0 (14) 


for all j, will be displayed and discussed later. 


Iterative Matrix Formulation 
The Sy algorithm (12), (14) with source (13) has the 


matrix description 
0-5 سپ(‎ )4+1[ fe sy (4) +q (15) 


derived from equations (11). 
If the equations are arranged so that the solution 
vector has the order 
92.1 
TN 
Phra 
92 1 
ý = "52 (16) 


የጄ 1.2 


۲,۵۶ 
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then 


ስ 2 (17) 


PIE 
mn ٤# 


The matrix (D-E) has the block form 


(18) 





N/ 2 


where each block has the lower bi-diagonal form 
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e a: = B. (19) 
© - ند‎ `d 
ፌ 


The index ZZ is that for the ordinate H... Each block is an 


R by R submatrix. 


The scattering matrix has the form 


(20) 





5: ia = Sgp. (21) 
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and 


O 
ው መመ —=— لیا‎ 


Observe that this matrix formulation is the upper left 
quadrant of blocks of the full-range Sy matrix formulation 
and only the elements 5j are different. 

Due to the block diagonal structure of (D-E], and 


the bi-diagonal form of each block, (p-E) ^! is simply 


EE OD) 





where 
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di 
e; 1 
d 

es ec ] 
E.g. a 

A A 

: -1 
k-1 at? = B. 
E 4 e è © * 99 © © © © © © © 

pth now EL Ei ۱ 

: : = 
ENS nenn. እ (24) 
38! 1 ር A 

: -2 7 
2 I EN Dus ° e T3 
ER .R-1 3 2 d. 
> ar dr dr ^ 


and 4 is the index for ur. Each block is of dimension R by 
K. (p-g) ^! can be applied in its general form to equation 
(15) to produce 

(x) 


(4+ 1) + (D-E)7?g . (25) 


- (D-E) ^! Sy 


Y 


In this equation 


(26) 


(D-E) 7's 
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Observe that this matrix is the upper-left quadrant of 
blocks of the full-range iteration matrix of the previous 
chapter. 


Each block has the form 





1 
TR 5 
° , | 
gom a 
2 ` 
BET. 7 
a? a? q? 
( í ( 
EE | 21 ٢ 
" 4 m i 
1 | 
k ^ QU al a l af 1 


Oo 
> 
' 
A 


Sala: 
Gs 


تما 


e Q 
a. aje. ع‎ 
' ( 
دحم‎ 


t 
دم‎ 





o 
1 


Calculation of the Norm 
The infinity-norm of (هدم)‎ 15 as previously defined, 


is easily calculated. The block row sum is 


E M -1 g^ = 
B; Spy +B; S E EBS = — (u, tw ٠۰۰۰7۶ ز‎ 8 5 
A BD, A BD, A BDy y 2 2 | ] 27 "7718 L 80) 
E 
a o! 
را سے ۔‎ Sgp (28) 


due to (9). Here 589 has the form of (21) but with unity 


I 4 
in each non-zero matrix element. Actually (28) is 25 times 





=> 


the matrix part of equation (27). Consequently 
^ e: R- 1 ] 7 er e ej s (29) 
2] _ O A LP A ee et 
ES | pos 5 2. 2፤ ) 

di 1 4 7 
where as before, 4 is an index for the minimum absolute 
ordinate value in the set {uj}. If one defines a f 

e. 
0 < cS RE j 
1 d: (30) 


equation (29) becomes 





4 R | 
x] O 1-8 RE 
D-E = — - 3 
|| (D-E)"*s]| 2 1 ) 8 | (31) 
ፈ 
Sufficiency conditions for ||(D-E)7?S|| < 1 are derived by 
requiring that 
Á 
c ] 
« ] 
ር 0 
er that 
2)) + Aot > 240° + I2 (u ;) - Nema (32) 
min min 


This is the same inequality which gave the sufficiency con- 
ditions in Chapter II for convergence of the full-range Sy 
algorithm to the exact solution y*. These conditions are 


the same for both problems and will not be restated here. 
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Exponential Spatial Transform 


The transform 


: (x] -$, Ix] e^* , (33) 


Y, j 


where a is a constant greater than zero, applied to the 
angularly discretized equations (8) gives 


Eo ( x] N/2 
መበ x T | — OX 
+ lot + Os) O; (x)= RAE + ۳۶ء‎ )34( 


۳۰٦,٠٠٠ ,N/2. Observe that the transform has added 
absorption to the problem in amount OM ¢ to produce a new 
"effective" total cross section, (ot + ak p) as well as 
changed the source. 

Spatially discretizing these equations in the same 


way as previously described for the unmodified equations 


gives 
(i+1) (itt) _ o? ÓN [4] [<] E 
AA > دج‎ ተክ رر ۷6ل‎ 7" Z (35) 
2 


K. R and í = 1,2,***,N/2. The algorithm in machine 


form is 


8 11 ፡/ [ረ-1]. ME n 3 - ox 
9b], j 74 j ‘AAT e Rt 1/2 y 3g) 
2 


. 


, 


where 


+۲ — Ld j (37a) 


and 





Bee 
(37b) 


Matrix Formulation 

The $-domain algorithm has precisely the same matrix 
form as the unmodified Sy equations. That is, equations 
(11) through (27) describe the -domain algorithm if ej and 
d; are replaced by e; and d; and if source and solution 
vectors are changed accordingly. 


The infinity norm of the transform domain iteration 


matrix is 


۱ qeu] 2 
٦۲۰۰٠ د و چاو‎ (i ل‎ “ደ ም ] ) (38) 


A 5 


where 4 is the same index as before. Defining 


٣١ 
0 > . اہ وج‎ (29) 


٨ 


The above equation becomes 


4 R 
_ he 
۶ "٣٣١ Td; É 2 ( Y. | (40) 


Properties of || (D-E)7?S|| 

The infinity norm of the transform domain iteration 
matrix is a monotone decreasing function of the acceleration 
parameter, a. For small o, the slope of ||l0-E)^!S|| with 


respect to q 1s steep but ‘decreases to near zero at q = a* 


where 
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在 
1 _ (imin Ag 


0 
5 1ዞ/ ) ዘፈዚ 


(41) 


which is the value of a when ei is zero. The norm value at 
ar is 


_ ho? (42) 
E شح‎ 


]])٤- 1۱| 
ቹ 7 
Ore A MAN 


From this point the norm slowly approaches zero as a 
approaches o. This behavior of |I (2-217 !s|| with ais 


verified by analysis of equation (40) in Appendix B. 


Optimum Value of a 
Appendix B reveals the behavior of | (2-817 !s]| 


for a> a* to be 


Moc sil- گے‎ E 8 [ለ5)፻ › ፻[ለ51" Y 
a min (4 + de] = 
where 
ao= a* + € (44) 
and 
0 > g£ < e° , (45) 
In the region where le < 1, gees | essentially takes 
its value at o*. Consequently very little norm reduction 


occurs in this region. This defines a practical “optimum" 
value for «a, that is, a = o*. Figure 3.1 illustrates this 


"optimum" value for a particular problem. Observe that the 
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norm is not further reduced for values of a > 2.00, at 

least to five significant figures. For this problem, 

0» 2 1.80924. Hence an Sy problem of the type discussed 
here has a practical "optimum" acceleration parameter which 
can quickly be calculated from equation (41). Use of this 
value will result in the largest practical norm reduction 


with accompanying improved iteration efficiency. 


Range of Utility of Transform Method 

From the previous discussion, it is obvious that 
De si] has reached a practical minimum Bu when €; 
is zero. Consequently the transform method can signifi- 
cantly decrease the norm of the unmodified 5, iteration 
matrix only when DEI S 1S non-negative, that is, when 


A < 2 uil min : (46) 


ot 


Otherwise e; 


peesenegative, which guarantees that e; is nega- 


tive for allo » 0. In this case, the norm has already 
reached a value which decreases very slowly with increased 
m. 

The transform method has the greatest latitude for 
norm reduction, if (46) is met, when ^ is chosen very small. 
Equation (41) reveals that a* is large when A is chosen 
small. A large «* decreases the norm of BUSES more 
than a small a*. Table V shows this property for a specific 


problem. 








> / 6፦ 


TABLE V 
a* FOR VARIOUS A 


Problem parameters: 


جح جرس جم حب 一 Ss SAE‏ سی 


|| [8-፻]!5| | at a" | slab width = 10.0 cm. 








of = 1.0 cml. 


0.69839 

có = 0.9999 cm 1. 
0.41903 

lu.) . = 0.23862 
0.20952 0 

|| (p-E)^!s|| = 0.9999. 


Table V shows that small A allows large values of 0 
to be used resulting in vast reductions in || (p-5E)~/S|| 
and improved convergence efficiency. The choice of a, how- 
ever is based primarily on the estimation of tolerable 
discretization error for a given problem. Consequently 
the transform method's utility depends upon factors other 
than those which most improve iterative efficiency. Once a 
^ is chosen, however, a* is easily determined. Its use 
provides practical maximum acceleration of convergence for 


that problem. 


Computer Experiment 

The half-range problem algorithm previously described 
in this chapter was run on an IBM 360/67 computer with the 
following problem parameters: 

Slab width = 10.0 cm 


R 


30 spatial intervals 


A 


0.355353 cm 


ot = 1.0 cm? 





>ሃ//፦ 


o = 0.9999 cmt 

Six ordinate Gauss quadrature set 

First collision source of unit strength 

Fractional iterative error less than ¢ = 1074. 

The improved convergence criterion, described in 
Chapter II, was used. A program listing and sample output 
for unmodified Sy are included in Appendix C. The transform 
method program and sample output are included in Appendix D. 

Table VI compares the number of iterations to con- 
vergence for various values of the acceleration parameter, 


O. 


TABLE VI 


TRANSFORM METHOD RESULTS 


Maximum Percent- 
Percent age Difference 
Decrease in Unmodified to 
Infinity] Iterations |in Number of Transform 
Norm to Converge Iterations Solutions 


0.99990 
0.95435 
0.89332 
0.80727 
0.69839 


0.69839 





Observe that a = 0 is the unmodified Sy result. The 


Symbol * signifies the optimum 0. For R = 30, the condition 





Lo 


— is met and the unmodified Sy iteration matrix is 


g 
non-negative. 


eh 
< 


Table VI reveals that the transform method accelerates 
convergence of the Sy algorithm for the half-range problem. 
The acceleration parameter, QA, decreases ej and increases 
d; so that each non-zero element of the iteration matrix 


ያ 
is decreased in value (see equations (26) and (27)). This 
produces a reduced norm and concomitant acceleration of 
convergence. 
Values of a greater than a* used here do not further 
reduce the norm or number of iterations to convergence. 
The dependence of convergence efficiency on the norm is not 


linear. Acceptable improvement of convergence efficiency 


is attained for values of ® well below at, 


Calculational Effort 

The calculational effort of a machine algorithm is 
based upon the total number of machine operations performed. 
Since additions and subtractions require only an infinitesi- 
mal amount of computer time compared to multiplications 
and divisions, only the latter are tabulated. Table VII 
lists the number of operations in each component part of the 
unmodified Sy, and transform method algorithms. Here n is 
equal to (R-N/2), á is the number of iterations to convergence 
of the unmodified algorithm and j is the number of itera- 


tions required for transform algorithm convergence. 
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TABLE VII 


COMPUTATIONAL EFFORT COMPARISON 


Unmodified Sy Transform Method 
Quantity Algorithm Algorithm 
Mesh sweep 
Scattering source £[N/2+2] 1 [N/2*2] 
Coefficients e, d 9(N/2) 13(N/2) 
Transform source 0 n 


Transform solution 0 n 


Total ፈ[28#8ታ2]*9 [2 /[28*##*2]*28| 2 [፡)] 





The transform method reduces computational work if, 


4 (3N-1 
DP E Lee C (47) 


For most problems R >> N so the transform method has utility 
even if only a few iterations are saved in its employment. 


53 : مر‎ : : 
Isaacson shows that Gaussian elimination requires 


+ n 


2o. 
n (n*-1) 2 (48) 


operations to solve a matrix system corresponding to the 
exact solution, y*. Comparison of the Sy iterative technique 
shows a significant economy for large systems, n large, if i 
is not too large. Even so, i would have to be of the order 
n? before the computational effort of the unmodified Sy 


algorithm reached that of Gauss elimination. 





E 


Error in Transform Method Solution 


The new convergence criterion 


NU g (4) 


J 
° 


08 9 ٦٣ 
«y [rediens (49) 


which guarantees that 


EEREN | 


Pp Pk 
v Et (50) 
Pk 
for all k, was used in the computer algorithm. Using the 


transform (33), inequality (50) implies that 


ype e 1 ver % 


፣ የከ oil — * E 
Vy, e 2 


or that 


— 118 (51) 


for all k, where yt) is the y-domain transform method final 
solution. Consequently the convergence criterion (49) and 
convergence properties of the method guarantee that the 
transform method wW-domain solution has fractional iteration 
errors less than €. This same criterion was met for the 
unmodified Sy, algorithm. The origin of the discrepancy in 


the p-domain solutions for transform method and unmodified 








-81- 


SN solutions, exhibited in column 5 of Table VI, is there- 
fore not attributable to iteration error. 

The origin of the discrepancy mentioned above is 
discretization error. The $-domain exact solution, $*, has 
much greater variation over the slab width than does the 
unmodified S, exact solution, y*. Figures 3.2 and 3.3 
demonstrate this for a particular direction. Since the 
spatial discretization of the derivative is accurate to 
order A* in both domains, it is a poorer approximation in 
the d-domain, where the solution has large curvature magni- 
tude, than it is in the unmodified Sy v-domain where the 
solution has a relatively small variation over the slab 
width. Consequently larger discretization errors occur in 
the transform method ¿-domain solution when large values of 
A are used. At 0r*, gi varied seven orders of magnitude 
over the slab width whereas the unmodified Sy v-domain 
solution varied by less than a factor of two. This accounts 
for large discretization errors in the y-domain transform 
method solution when Q is near a*. Table VI reveals, how- 
ever, that as à is decreased from a*, the discretization 
error diminishes so that the y-domain solutions for both 
algorithms agree to within less than 0.4% at a = 0.2. For 
this A, the d-domain solution varied less than a factor of 
five over the slab width, producing smaller discretization 
error than did larger 0. Here convergence was accelerated 
to 78.5$ of the maximum acceleration obtained at a*, The 


value of a = 0.2 was a good practical choice for this 
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SAMPLE UNMODIFIED Sy p -DOMAIN SOLUTION 
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FIGURE 3.3 


SAMPLE TRANSFORM METHOD 中 -DOMAIN SOLUTION 
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problem. For other problems, the user must make a choice 
of a based on practical experience and a guess at the 


expected solution. 


Influence of A 

Previous discussion connected with Table V revealed 
that, in theory, a smaller A allows the choice of a larger 
a* with accompanying larger decrease in ||(D-E)~/S|| ana 
number of iterations to convergence. In practice, however, 
such large values of ® are not useful because of the 
extremely large resulting discretization errors in the 
$-domain solution. Another consideration also limits the 
practical usefulness of large 9. For thick slabs, the trans- 
formed source, which has a factor وت‎ can become extremely 
small. Since all computers have a limit on the smallest 
number it can handle, E for the IBM 360/67, it is easy 
to exceed this limit for large o. 

As a consequence of the above factors, it is best to 
use values of a much smaller than a*, for a problem with any 
A for which the transform method is applicable. Computer 
experiments conducted by this investigator have verified 
that the choice of & = 0.2 produces approximately the same 
acceleration of convergence and accuracy listed in Table VI 
for the problem described previously but with the much 


smaller values of A resulting from R = 50 and R = 100. 





= 
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Generalization of Transform Method 

The transform method immediately generalizes for the 
half-range problem treated here to include heterogeneous 
slabs. In this case the cross sections are functions of 
position. The relatively simple expressions for 0 Si] 
and ||(?-E)^'S|| shown in this chapter do not occur but the 
matrix block form prevails and these norms can be calculated. 

Transforms such as equations (3) and (4) place space 
dependent absorption into the problem. When using trans- 
forms of this type, care must be exercised to restrict the 
range of the arguments in order that finite non-negative 
absorption is added to the problem by the transform. Other- 
wise reduced or negative "effective" total cross sections 
will result and convergence will be decelerated in the 
transform domain. 

It is possible that problems in 4,u and higher- 
dimensional geometries can be formulated in the mathemati- 
cal context displayed here for the simplest geometry. If 
transforms can be found which place an acceleration parameter 
in the proper place in the iteration matrix to reduce its 
norm, vast savings of computational effort can be expected 
for these problems. The transform must be one which leaves 
the angular discretized S, equations invariant. To be of 
practical value, such a transform must not only properly 
position an acceleration parameter in the iteration matrix 


but also must not alter the source so drastically that large 
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variations occur in the d-domain solution with the 
accompanying unacceptable discretization error. 

Possible transform candidates are some form of the 
integrating factor for the angular discretized equations 
describing the problem after neglecting the scattering gain 
term. Transforms are not restricted to this class, however, 
and will prove useful if they meet the criterion described 
above. 

It is doubtful that a transform can be found which 
is useful for all problems in a particular geometry. 
Experience of this investigator indicates that it is more 
likely that each problem in the geometry has a "best" trans- 
form which can be prescribed based upon the source, boundary 
conditions, and the prescriber's estimate of the form of 


the solution. 


Conclusions 

1) The transform (33) applied to the half-range slab 
problem accelerates convergence of the Sy algorithm for 
predominantly scattering media using the improved convergence 
۲۶۱۹۴0۰ 

2) The transform method is ED only for prob- 

ን 2 un 

lems in which ^ is chosen so that A mo 

3) Significant convergence acceleration occurs for 


values of a much smaller than a*. These smaller values of a 


produce solutions with acceptable discretization errors. 
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4) A practical choice of 9 for a particular problem 
must be based upon experience and the user's estimate of 


the expected solution characteristics. 





CHAPTER IV 


DISCRETIZATION ERROR IMPROVEMENT USING 
SPATIAL TRANSFORMS 
Discretization error or truncation error, as it is 
frequently called, results from discretizing an equation in 
the continuous domain and comparing the discretized equation 
to the analytic equation. For the full-range Sy approxi- 
mation of Chapter II, this process is illustrated by 


approximating the set of ordinary differential equations 


d .[ x) £ 4 N 
EO 
“መ” ہہ ٭‎ Len ll + q;lx) a) 
j = 1,2,--:,N, which are valid in the spatial continuum 


ሀ. < x < L, by the spatially discretized set 


Pam + teen ji) 5 \ | ም a , 
Se) °° s — ا‎ Ga een i to) 


17 


which uses a centered difference approximation for the 


derivative and a simple average for y , at each mid- 
Xp, 1 *Xp ert 


point X12 A سے‎ 
Expanding Vi, 7) and TE in Taylor series about 


the point VEU و‎ truncating appropriately, and forming 


the proper sums produces 


لو ا ان 


V; Xp, 1 yo) * V" (X,,1551 0la*) — (3)‏ . فل 
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where X, < X.,1/2 < Xp +] (4) 


and 


Vi (Xoo) th; (Xp) 


ጠጠ O || ۷ اس‎ "Iv 
where X, « ہرمک < ما‎ : (6) 


Substituting equations (3) and (5) for the appropri- 
ately indexed quantities in equation (2) gives 
۸۸1 21 , vr, پڑت‎ . M ር 2 "1 : የ1". 1 gn oie N 
2 2 4 
(7) 
Now evaluating ۷ [x] at X +1/2 in equations (1) and sub- 
tracting the result from equation (7) glves the discretiza- 


tlon error 


N 
E 2 im f 2 n= _ 2 " ws 


If it is assumed that the solution y LXI to equation 
(1) is sufficiently smooth that its third derivative exists 


and is bounded, that is, 
y'"(x) M (9) 
J 
for all j,x, then 
| « MO(A*) + pO [A2] (10) 


where 


eos N کے , ا‎ 
0 = nax eta) - RA) (11) 
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The discretization error is clearly of order c for this 
method. That is, T approaches zero as Ke approaches zero. 

A standard technique for diminishing the discretiza- 
tion error is to choose A small enough, depending upon one's 
estimate of the variation of the solution, that the dis- 


cretization error is tolerably small. 


Transform Method 

The above technique of diminishing A to improve the 
discretization error has a practical limit. As A decreases, 
the number of spatial mesh points increases as does the 
computational cost. Consequently a choice of A must be 
based on the tradeoff between accuracy and cost. 

Another method of decreasing the discretization 
error is suggested by (10) and (11). If a transform can be 
found which renders equations (1) invariant in a transform 
domain and diminishes the second derivative of the solution 
without increasing M, the discretization error will be 
reduced. That is, if the transform domain solution of 
equation (1) has less curvature magnitude over the slab 
than the solution ALS of equation (1), the transform 
method will improve the discretization error. This method 
Will be demonstrated in the next section for a simplified 


Sy problem which has an analytic solution. 


Simplified Problem 


The equation 


a + oty (x) = o plx} + qix) (12) 
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with boundary condition 


VIG T0 

and source 
- tx 

q(x) = e 

has a solution 
“ر‎ (0% y 0-1) 1. 
በከመ. 2 477 di ^ 
(c -ugO | ا‎ 

where 

ot = g-o 


and u, is a positive real constant. 


(13) 


(14) 


(15) 


(16) 


An Sy-type iterative solution of the same problem is 


formulated by approximating the derivative by a central dif- 


ference and approximating equation (12) at a spatial mid- 


point as before to give 


vit 
PESE e (itl) k+ 
Ub , 1 = J Vk (nr 
with 
V 2-9 - 
and 
_ 


vat!‏ و 


* 05,12 


(17) 


(18) 


(19) 








دوو-۔ 


Here 
Zu, As 
e = 一 — (20a) 
and 
Zu, * Ag 
Q = u De (20b) 


Observe that this algorithm is precisely the same form as 
the Sy algorithm along a Iu; ray with the exception of the 


quadrature approximation of the scattering gain term. 


The Exponential Transform 


The transform 
wix) = olx) e ^, (21) 


with a » 0, can be applied to equation (12) to give the 


transform domain equation 


uy, A * (o*- ou, 19 x] = o? $ (x) j qix) e” : (22) 


Discretizing spatially, as before, gives 


ES 
(61) — (امغا ع‎ 1) a( he] * dk اه‎ 
9 + 1 d 9k = ane || لا‎ En; — (23) 


2 
where 
?u,-Ac «ay, 
e = ES (24a) 
and 


29. ተለዓለ... 


d = ESSE 2 (24b) 
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Observe that this transform places the parameter % into 

such a position that B is increased and = is increased. As 
shown in Chapter II, this increases the norm of the iteration 
matrix and decelerates convergence of the algorithm. Con- 
sequently the method is restricted to those problems for 
predominantly absorbing media in which the iterative 


efficiency is not so sensitive to the norm of the iteration 


matrix. For these problems, the conventional convergence 


criterion 
إن‎ 7 Tie 
وو سب‎ a - < € (25) 
ሃኔ 


can be used with negligible loss in iterative precision. 
The exact -domain solution of equation (22) with 


source 


ax - «t a-g* 
o NONE 0 SUUS s E e (26) 
and boundary condition 
p[o) = 1.0 (27) 
EX =) (o-u o*-1) „le -au,] 
is ከ. . -.ቭሽ !፦ <“ ፡፡.. . - (28) 
CS lou at) 


Sample Problem 
A problem with the parameters listed below was 


treated by the methods described above. 





=94- 


Slab width = 10.0 cm. 


R = 30 spatial intervals, 


gt EIE era 
o = 0.2 ne 
N 0.9. 
e = 1074, 


The exact solution of equation (12) for these 


parameters is shown in Figure 4.1 and is 


a ~ = 400 . (29) 


The exact $-domain solutions for various values of o are 
shown in Figures 4.2 through 4.6. For values of a = 0.5 
through 0.9, the exact d-domain solution has much less 
curvature magnitude over the slab width than does the exact 
v-domain solution. Consequently it is expected that the 
transform domain iterative algorithm (23) solution will 
possess smaller local discretization errors than the 
unmodified Sy algorithm (17). 

Algorithms (17) and (23) were run on an IBM 360/67 
computer to verify the expected discretization error 
improvement using the transform method. Table VIII sum- 
marizes the results. 


The fractional errors were calculated from 


k j uf (30) 
exact 
R 
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FIGURE 4,1 


EXACT SOLUTION FOR UNMODIFIED Sy 





8 
& 
F Y (x) = N 
5 
1 
፳ 
8 
. 00 2.00 4.00 6.00 8.00 10.00 
X, CM. 
FIGURE 4.2 


EXACT SOLUTION FOR a = 0.5 


s 中 (X) = > - Mop S9 
š, 
i 

8 

2 

P 2.00 4.00 6.00 8.00 10.00 


۳ھ 





NORMALIZED SOLUTION 


NORHALIZED SOLUTION 


سنا 


.00 


— 


69 & 


e 


.00 


-96- 


FIGURE 4.3 


EXACT SOLUTION FOR a = 0.7. 


٠۔ا ,رم‎ ٦٣٤ _ 10e ጓዳ 
2.00 4.00 6.00 8.00 10.00 
X. CM. 
FIGURE 4.4 


EXACT SOLUTION FOR a = 0.9 


ES _ Ex 


Di x 11 10e 


2.00 4.00 6.00 8.00 10.00 
X, CM. 
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.00 


69 


80 
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«20 


.00 


ey m 


FIGURE 4.95 


EXACT SOLUTION FOR a = 1.1 


$ lx) 5 1001 7 .م70‎ 


2.00 4.00 6.00 8.00 10.00 
X, CM. 
FIGURE 4.6 


EXACT SOLUTION FOR a = 1.3 


lx) » 11e: 4! * - 104:?X 


.>-፦ጅኤጅኤ ልመ 


2.00 4.00 6.00 8.00 10.00 
X, CM. 
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for the unmodified Sy, algorithm and 


exact (i) 
>, E | (31) 
exact 
5 


for the transform domain algorithm. The errors compared, 
however, are for values of the unmodified Sy ~-domain solu- 
tion and the transform method y-domain solution. These 
fractional errors are attributed primarily to discretization 
error because the unit roundoff error for the double pre- 
cision mode of the computer is 10714 and the iterative 

error should be less than 10”? due to the convergence 
criterion (25) used. Certainly fractional errors greater 
than 1073, as most unmodified Sy errors are, can be attri- 


buted to discretization error. 


Discussion of Results 

The transform method produced steady discretization 
error improvement with a for all space points up to a = 0,9. 
This error improvement is due to the general decrease in 
curvature magnitude of the exact d-domain solution compared 
to the exact yU-domain solution over the slab width. The 
transform method can be used to diminish p in equation (ll), 
thereby reducing the discretization error. 

The value 9 - l.l produced error improvement over the 
unmodified Sy Solution but less improvement than did a = 0.9. 
Values of A > 1.1 produced larger discretization errors than 


did unmodified Sy in some regions of the slab. Observe that 





0/0 — 


the solution $(x], Figure 4.4, has the least curvature 
magnitude at % = 0.9. Also at a > o, = 0.9, equation 
(28) reveals that $(x) has a growing exponential in the 
second term Ens dominates the solution. This suggests 
that predicting an optimum value for Q requires knowledge 
of the exact solution. Obviously this is not possible in 
more complicated practical problems but some guidance is 


available from the simplified problem. 


Strategy for Determining Optimum oa 

A practical value for a may be obtained from the fol- 
lowing argument. The method is only useful for predominantly 
absorbing media. These problems are characterized by solu- 
tions which "follow the source." That is, if the source 
distribution is a decaying exponential, the exact solution 
will have the form of a decaying exponential. Consequently 
a transform which changes the source to nearly a constant 
in the $-domain would probably produce a more slowly varying 
$-domain solution with accompanying smaller discretization 
error. For the problem treated here, that would be 


1 
q(x)e = ue ee const, 


when 


This is very close to the value o/u, which was found nearly 


optimum by knowledge of the exact solution. In the absence 





= 101 - 


of knowledge of $(x], the above strategy based on source 

shape could be useful. Of course if partial insight into 
the expected solution is available from analytic solutions 
to similar problems or if the user's intuition is sharp, a 
value of 0 based on such considerations could prove fruit- 


tul. 


Application to Sy 
The Sy algorithm for the full-range slab problem, 


derived in Chapter II, is repeated here 


4 
O رس‎ 


EE. -í (4 | 
Y Y کہ‎ vit f q : 


+ 
k+l,j d; 1 


QS uS uen uy d je I o E rg 
b, j d; "Rej - 05V br1/2,n ፡ በ رر‎ 
for Aj < 0 
where 
? |u ,| -Ac* 
5 ከመ — 
and 
2|u .|+A0% 
HE m ee ves (32d) 
1 2A 
The transform 
۱ اد‎ PAX 
ሃ/!።) 0 (21) 


with 4 > Q applied to equation (1) gives 





do . (x) ፅ አ 
7 — + (ot an 19 00) = TI በቀ !።] ۴ ajlde” (33a) 


for H; > 0 and 


Á 


dd; (x) 
ul + total, ie سو« (ج)‎ 


for Hy < 0. Observe that the transform adds absorption to 


N ax 
mum ode (33b) 


the y < 0 equations but diminishes absorption in the 


u; > 0 equations. 
Spatial discretization is conducted in the manner 
previously described and coefficients of the fluxes are col- 


lected to give the Su algorithm, 


Q X 
e R+1/2 (34a) 


[4] 
ne +1/2,n 8+፲/፻,/ 


1) 1 Jos N 
m 
n= ] 


ES €; (4 
Peel, j Tk, d; 


for u; > 0 and 


j 
EET *; (4*1) ,1 Jo? VE OX 1/2 
"hei وت‎ 8712 Donee 12 n 9112, (34b) 


for ur < 0 where 


"|0 ار‎ 
IE — for ولا‎ < 0 (35a) 
alu; 


$ 
= ~ for Hu; < 0 (35b) 
and 

a lu - | 

al E — for u> 0 (36a) 
d: = 

j AA 

dj + =p for uj < 0 (36b) 
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Observe that QU increases 7 and m for uj Re 
Consequently the part of DEIS which corresponds to the 
7 > 0 equations has a positioned such that each non-zero 
matrix element is increased. Conversely, Q diminishes í 
and = for ug < 0. This results in a decrease in the non- 
zero E elements of the iteration matrix corresponding 
to the uj < 0 equations. These properties destroy the 
simple block collapsing technique, demonstrated in Chapter 
II, used to calculate the norm of the transform domain 
iteration matrix. Some portions of the maximum row sum are 
diminished and some are increased by 0. Fortunately it will 
not be necessary to calculate E) CN for the problem 
treated here since only predominantly absorbing media are 
considered. Consequently the norm will not change signifi- 


cantly with nor will the number of iterations to con- 


vergence. 


Sample Problem 


Equations (34), subject to the boundary conditions 


ቀ 


102 for h, > 0 (37a) 
La J 


and 
Am. tor cQ (37b 
unm; u; ) 


and source 


q [x] = 0 for all j, (38) 
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were solved on an IBM 360/67 computer. The pointwise con- 


vergence criterion 


(i+1) (x) 


HU = 


4 
% 


< € (39 ) 


was used. The problem parameters are: 
Slab width = 10.0 cm. 
R = 30 spatial intervals 
o* = 1.0 ml. 


o® = 0.2 mt. 


Quadrature set = six point Gauss-Legendre 

a = Various values from 0.1 to 1.0. 
The resulting y-domain solution will be compared to two 
unmodified Sy solutions of the same problem, described by 
equations (32), but with different spatial mesh sizes. 6 
used R = 600, the other R = 30. Appendix E contains the 


programs and sample output. 


Unmodi fied Sy Results 

Comparison of the unmodified S, results for R = 600 
and R = 30 shows that all R = 600 values are either equal to 
Or greater than (at least to five significant figures) the 
R = 30 solution values. Since the R = 600 solution has less 
discretization error than the R = 30 solution, the above 
fact simplifies the analysis of the effect of the transform 
method on improving discretization error. If the transform 
method Y-domain solution values are between the unmodified 


Sy R= 30 and R = 600 values, the method has improved the 
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discretization error. The above criterion is used to judge 
the effectiveness of the transform method in improving 


discretization error. 


Transform Method Results 

The data show that the transform method decreases 
the discretization error for the problem. Discretization 
error improvement is not pointwise uniform but, in general, 
transform method solution values approach the R = 600 solu- 
tion values as o is increased. The bulk solution differ- 
ences, displayed and defined in Table IX, diminish for each 


u, component of the solution as QA increases. The total 


1 
bulk difference also diminishes as O increases. This means 
that some integral measure of the discretization error 
diminishes as a increases through a = 1.0. A closer look 

at the pointwise data, however, reveals that at a = 1.0 

many transform method Y-domain solution values have overshot 
the unmodified Sy R = 600 solution values. 

Table X displays a more detailed analysis of the 
pointwise discretization error improvement achieved by the 
transform method. It shows that discretization error is 
decreased at each point for values of & up to 0.5. For 


a > 0.5, discretization error is generally improved but 


some values overshoot the R = 600 values. As à grows larger, 


more values overshoot the R 600 values and the analysis 
loses its basis for comparison. All methods required the 


Same number of iterations to converge; 10. Hence 
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discretization error improvement was achieved at the expense 


of only R-N extra computations. 


Conclusion 

The transform method can be used to decrease the 
discretization error in Sy problems. This method is effec- 
tive when the transform domain solution has less curvature 
magnitude over the slab width than the unmodified Są solu- 
tion. 

The transform used for this purpose must be tailored 
to each individual problem. The choice of a successful 
transform for a specific problem is based upon intuition 
developed for that problem by experience, knowledge of 
analytic solutions to Similar problems, or guess. For 
problems in which solutions tend to "follow the source," a 
transform may be devised which converts the source to a near 
constant in the transform domain. A possible transform 
candidate is some variation of the integrating factor for 


the transport equation neglecting the scattering gain term. 





CHAPTER V 
ROUNDOFF ERROR ANALYSIS FOR ITERATIVE METHODS 


Background 

Previous chapters have discussed the concepts of 
iteration and discretization errors in the Sy method. Each 
effect was treated with regard to its role in the con- 
vergence of the algorithm to some exact solution. The 
approximate solution Se! was shown to approach the exact 
solution y jx) of the angular segmented equations as 4 
approaches œ% and A approaches zero. This concept of con- 
vergence ignores the errors that occur in any machine 
Computation due to the rounding procedure employed in 
floating point arithmetic. In order to be effective, how- 
ever, an algorithm must remain immune to the accumulation of 
roundoff errors. This immunity is called numerical sta- 
bility. An otherwise convergent algorithm may be driven 
divergent if it does not possess this property. In order 
to analyze the property of numerical stability, precise 


definitions must be made. 


Well-Posed Computation 
A numerical method, which produces a solution 
approximating the exact solution of the problem, is said to 


be a well-posed computation if the solution is insensitive 
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to small changes in data or to roundoff errors. Here data 
refers to the elements of matrices involved and the initial 
or boundary conditions and sources. Isaacson?? defines a 
computing problem as an algorithm or "a set of rules 
Specifying the order and kind of arithmetic operations 
(i.e. rounding rules) to be used on specified data." He 
also defines a computing problem to be well-posed if and 
only if 
i) a solution exists 
ii) the solution is unique, that is, when performed 

several times, with the same data, identical results are 
obtained 

iii) the solution depends Lipschitz continuously on 
the data with a constant that is not too large. The last 
condition requires that "small" changes in the data must 
result in only "small" changes in the solution. A well- 
posed computation possesses the property of numerical 
stability. A necessary condition that a finite difference 


scheme be convergent is that it possess numerical stability. 


Condition Number 
Roundoff errors will in general be introduced in the 
course of carrying out the arithmetic required to solve 


the system of equations 
Ax = í . (1) 


But if the algorithm is well-posed, these errors can be kept 


within reasonable bounds. The matrix Á is said to be 
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"well-conditioned" or "ill-conditioned" if the computation 
is or is not, respectively, well-posed. The condition 
number, v{A), for A serves as a measure of ill-conditioning. 


The condition number is defined by 


v{A) = | 0۶۰۰۰ (2) 








Some important properties of the condition number are: 


4) ህ[ላ] > 1 with equality only if A is orthogonal. 


ii) v{A) = v{A7!). 


iii) v(«A] ። ህ [ላ). 
Suppose that the data A and 4 of equation (1) have 
small perturbations or uncertainties $4 and 64 respectively 
due to roundoff error. The computer is solving the slightly 


perturbed system 


(A + SA)(X + éX) = í + $í (3) 
instead of system (1). Isaacson? shows that if 
1! ዚዜ 1. ll, (4) 


the resulting relative uncertainty in the solution, X, is 


| 16: | 5 V P) 
HxH T S 0 ጠህ '6ባ 
||ለ|| 


Observe that if v(A) is relatively small, A is well- 
conditioned and "small" relative uncertainties in the data 
produce "small" relative uncertainties in the solution. 


Conversely if v(A) is very large, A is ill-conditioned and 
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"Small" relative uncertainties in the data can destroy the 
solution. 


64 points out a popular misconception that 


Forsythe 
the smallness of det{A}) causes the ill-condition of A. 
That is, the condition number is not a measure of how nearly 
A is singular, hence difficult to invert. He stresses that 
the size of v(À) is a far more important criterion of the 
"badness" of computational system (1) than either the 
smallness of det(A) or the largeness of the order n of the 
system. 


A theorem from Isaacson, which will prove useful 


later, is stated here. 


Theorem 2 
If the matrix A has ||A|| < 1, then (I + A) is non- 


singular and 


1 


ር SRE ee \ 
+1 2 


= ] 
2 < TTA ` 


Application to Iterative Methods 


33 referring to iterative methods for 


Isaacson, 
solving systems of equations (1), states "One of the intrin- 
sic advantages of such methods is the fact that errors, due 
to roundoff or even blunders, may be damped out as the 
iterative procedure continues". The following analysis 


provides an insight into the mechanism of this intrinsic 


advantage. 
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ein terative Sy algorithm, described in Chapter IT, 


has a matrix representation 


+08 ا‎ = sy) ta (7) 


At each iterative step, the data D, E, S, q, and y C are 
known and the machine is expected to invert (D-E) exactly 
and multiply the result onto the right-hand-side of (7), 


that is, solve for 


EET) (x) l 


ህ ር ር ከህ per s (8) 


precisely. The computer, however, has a limited storage 
capacity for each digit and must round off the input data 

to the last significant digit of its capacity. Consequently 
the uncertainty in the stationary data is ôD, ôE, ôS, and ôq 
which are of the order of the unit roundoff error of the 
machine (10714 for the IBM 360/67 in double precision mode). 
These uncertainties are propagated through the arithmetic 
process of inverting (D-E) and matrix multiplying the 

result onto appropriate right hand side vectors. Conse- 
quently even the input data اس‎ has its uncertainty, sy), 


carried with it in equation (7). Therefore the exact system 


to be solved at the i+ prt iterate is not (7) but 


011 (9) 


(D-E) y V. 


where the previously computed solution, ve, 15 


oyli | )10(‏ + لكان . لكان 





د٢‏ ته 


But now the machine must solve the perturbed system 


(8+85-5-65) (yp isp 4T) = (5+65) ray, +(q+6q] (11) 
where 59.“ is the roundoff error resulting from machine 


[4]. 


handling of the data Ya 


It is the same order of magni- 
tude as ôD, E, ôS, and ôq. These uncertainties in the input 


(<+7) in the solution. USing 


data produce an uncertainty OY 
equation (ll), it is possible to derive an upper bound on 


the norm of the relative uncertainty in the solution 


(4+1) 
ehr T at each iterative step. 


Flrst Iteration 


Suppose the initial guess for algorithm (7) is 
yo! ig (12) 


with no machine error. That is, the first machine iteration 


is on -m The first iterate is 


5-8)! = 0 (13) 


but the computer solves the perturbed system 


E easy!) = o 6 8g. (14) 


Factoring out (D-E) from the left-hand-side and multiplying 


by (D-E) l gives 


Md = (D-E) ^! (q*6q) . (15) 


1(5-5) ”! (85-65) 0 LE 
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Recognizing that 


- ] 
p! | = (D-E) q (16) 


from (13), cancelling this equality in (15) and collecting 
(1) 


terms in ôy on the left-hand-side produces 


=] 


[፲+(8-8)  !(59-6)| 69! » (5-5) ”! (68-65)ሀ!!፤፥ (5-5) '64 . (17) 


Hence 


aya [1+ (D-E) 7! (6D-6E) | (D-E) (8E-6p) y!!)4(p-E) 7489. (18) 


But due to Theorem 2, 


21 -1 E 1 
A 2 -፦፦------ፐ--፦ (19) 
1- | | (D-E) ' (65-56፻) | | 
subject to the restriction that 
EEE (20) 
Press si | 
Hence 
||[69'"!|| NE ۱ + [EU 
1- | | (D-E) 7! (6D-6E) | | ۱ 
Now divide by a | | and multiply by | = to obtain 
] 6D-6E Š ú 
lis j|. ሠ ከ.“ “ሆመ (21) 
14. ፤-|[| ው)! (6-88) || || ወ-8) || || -8) || !||1!? || 


where v is the condition number of the matrix (D-E), that 


is, 
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v = | | (9-5) || || (5-5)”!|| . (22) 
Bue from (13), 
_ (1) 
٣٣٣٠٠٦٢ ۱٢ 
Or 
. 1 (23) 
| [D-E] Jella Il ۱۱١١ 








From restriction (20) and the fact that 








|| (5-5) ”፤ 65-65) || « || (5-5) !|| “| |ፅው-65|| < 1 , 


|6D-6E|| hence 








I- || «»-E) ^ t6p-8E) || > 1-11 (0-8) 7? | 


] 1 
7 : 5; G 
1- | | (D-E) ` (8D-6E) | | |- || (0-5) | | SD-6E| | 








Using (23) and (24), inequality (21) becomes 


P] 
Eh. < 


E 


ህ | | SD-6E | | 84 


Saas ar EIE 
| ٢ 8 NS ||51| (25) 
D-E ۱ 


ہ 


Second Iteration 


The second iteration treats the problem 


(5-ጅ)ህ') = sy.’ + 9 (26) 


but the machine solves the perturbed system 


8121) = (8+65) (ሀ. oy!) 


— 


Bi 


(D+6D-E-6E) (v }+(g+6q). (27) 
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Following the same procedure used during the first itera- 











tion, 
el | 116591. ||*|| ۱۱''هٰ۱2|‎ , ۱۱۰۰۱۱ ۱۱ Eu (20) 
010068 Ces en Hoze jg T Toe ሞም |ውጆ|| 11"[| . 
[eal || 
A bound on the quantity —————— can be derived from (26). 
(2) 
B... ااا‎ 
-1. (1) EI 
y! د‎ (D-E) Sy, * (5-5) و‎ 
01 
=] 
] = 2) -] 
uU! - [o-m^s] [sy ?! - om "el. 
Consequently 


1" |. [ው '|1ህ""'ቨ » HE om HL a 


OL 


(7) 
lis, ال‎ || [mms] | ] + 116-8) 11-11411 (29) 
|| | lly || : 


A bound on ገ || on the R.H.S. of (29) can be derived by 
using the unperturbed system 
2 E I 5 
y! | = (D-E) Isy! | + (D-E) r : (30) 


From equation (13) 


yll = (D-E) lq (31) 





SS 
which substituted into (30) gives 
[I + C = ya 
or 
"1 = |] + (p-Ey !s| ١ | (32) 


Using Theorem 2, 


)2 
۳۷یپ 11 


۸۰۰۰۰٦ 
subject to the restriction that 
(8دم) زا‎ ۱ 2-7 NE )33( 


But recall that (33) is the condition that the Sj algorithm 
must meet to guarantee convergence so it is met a priori. 


Consequently 


] | 


የጠ s 5 “ገ But from (23), 
፪) || 110555) S| [+ 


ሀ!!! || 





so the above inequality becomes 


: 
| A 
هم‎ 


] > | |5-ጄ| | | (34) 
||ህ!"?|[|. ||2||[ EU 
y 9||”[፣-|| ወ-ሄ) '51|] 





Using (34) on the R.H.S. of (29) produces 


(1) š _ 
IL < lloms] I] pakez], c 
: | 


= 
|| ]-|| (5-5) 41|| 








= 








(1) 
Also using (34) for the Factor | Se [| in (28) provides 
the bound 2 | 
||69!!]|| Ilo || 
ሼር መጨ“. ር : (36) 
ll ll ||ፈ||[፣-|| መ-5) !51| |] 





Now applying (34), (35), and (36) to (28) produces 


se دا‎ EN: EE E EI H), 
1፡1. LS [፣-|| 5-5) !5| |]|||»-ጾ1|. || መ-5) Isl 


| Yes هاا‎ M الا‎ (37) 
- In 


where 


o - |[o-m^!spp- qp [o-m7!s] || (38) 








is the condition number of the iteration matrix. 


The General Iterate 


th iteration has been completed. That 


Suppose the 4 

is, pi) = yl) + sy!) has b ted. A bound h 
ae =ህ oy as been computed. ound on the 
relative uncertainty in the i+14* iterate solution can be 


derived from recognizing that instead of solving the system 


g (39)‏ لوو لی 


(D-E) y 
the machine solves the perturbed system 


(D+8D-E-6E) (y!4* Magy 447) - (segs) (yl4) sey!“))+(q+8q) . (40) 


Proceeding inthe same fashion as in the previous section, 
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[<] 
| دل‎ 11sE-6D|| , 5ة||‎ llego ll أا .ا ا‎ 
an eenn መከ Meng መሸ eng 
115311 
جار سر رر ڈو‎ (41) 
Ilo-zl]- | ly 7} 


From equation (39), 


70 — (D-E) sy (i + (5-5) '9 











Or 
pie) + [ሠ-5)!] [ptn 9 [D-E) lq] 
Consequently 
Iri ET | 5-5) '||-||4|| 
05 ፈቱ | ies || ው-፡) s| ü ] + +89 ” (42) 
The factor = |^ the R.H.S. of (42) can be 


|ህ 


bounded by recalling that the unperturbed system obeys 


; ፈ 
oe = } A" (p-g) ^! q (43) 
m=0 


where for convenience of notation 
A = (D-E)~/s , (44) 


Consequently 


፪። 6-8)7 04፦- E ty 





But if condition (33) is met, 


So 
۶سس۰‎ 4 - a yT mae 
= Mean - AS ۴۹ 
9 1 ا تت‎ 
٩11 ہی‎ ۰٠7۶٦ 
Hence 
q 7 + 70 
Consequently 








Hall < | [D-E | 7 ٦ wen 
E t- [ae] 


due to Theorem 2 if 
|14፤| ፡ | 
which is certainly true if condition (33) is met. Since 


Moss e MEE 


7+ ۶ ۶۹۹۹ ۹ ۹ ۶۵۶۴۱ 





ን 2 5- 


and 











, (46) 
AS TN ee 
AS 
Using (46), inequality (45) becomes 
7 ۶۸ء7‎ : (4+1) 
Ilall < هاا‎ 98 20 | | 
< E 
hl 
or in terms of the iteration matrix norm 
a. ] oe] ]-Er+] | (0-2) لاد‎ 
l Ilal > [1-1] 0-5) 7151147] (47) 
Applying (47) to the R.H.S. of (42) gives 
ዘበ 1. .ፒ ቦዩ. len SNN سم‎ "stri? 
Ip" || 6-5)” '5|| el diesen ٢ (48) 


Now applying (47) and (48) to (41) produces 


(4*1) " 
| |» || A ህ በ 1 2 11 177 7 | 
A ee ا‎ no cL ue ونت ال‎ 


Doc n 5... ||9-8||. ፤]-|| (5-8) ”፤5|!”ቫ| | [9-8|| || (9-5) !5|| 
| ነ 
USS [liom sf lel ° [rd om sl] 115911 

TaT MU UE (49) 


This expression is valid for 4 = 2,3,-:: 


Analysis of Error Bound 


Observe that in expression (49), the only terms which 
are affected by the number of iterations are || psg) | 


and ] 


A __. As 4 increases, the former magnifies 
1-|| )(9-2(-15| 1 





ከ دا‎ 


the bound slightly and the latter diminishes the bound. The 
latter effect is predominant. Therefore the bound becomes 
tighter as iteration proceeds and roundoff error propagation 
is diminished. This effect is not strong enough to force 


the relative uncertainty to vanish, however, as 4 approaches 


oo 


Certain terms and factors in (49) can be neglected 
with respect to more dominant quantities in order to obtain 


a qualitative view of the bound. That is, 























|1'""!! | : EE esl] «م‎ BINH a (50) 
Ber x دہ ہج کہ ےہ‎ š 


. 0 ¡- »+]]5D-5E 
D-E 


| |p-E| | ||o-E]| || ወ-5)!5|| |۱ Hall 


contains the dominant quantities in (49). 

The condition of the matrix (D-E) has a dominant role 
in the roundoff error bound. If v is very large, the rela- 
tive uncertainty in the iterate solution will be significant. 
If v is large enough that the factor | | is close 
to unity, the relative uncertainty Can be extremely large. 
For this to occur, v would have to be of the order of 1014 
for the double precision mode, extremely ill-conditioned. 
This event is unlikely. Since v, which is always greater 
than unity, appears to the second power in (50), it is the 
dominant quantity in the bound. 

The condition number of the iteration matrix, p, also 
plays an important part in the roundoff error bound. Its 
effect is subordinate to that of v since p appears only to 


the first power. But if (D-E)7?S and (D-E) are both 





5 


Lo could cause significant 


ill-conditioned, the quantity v 
roundoff error propagation in the solution. 

The error bounds derived here display another 
important effect. Expressions (25), (37), and (49) each 
contain the uncertainty ||6D-$E|| in important positions. 
This can cause the exact cancellation of roundoff errors 


during the iterative process. That is, ||$D-6E|| vanishes. 


This gives vivid corroboration to Isaacson's observation 


quoted previously. The presence of ||sSD-SE|| in the 
: ህ : 
dominating factor Kz Tép-sg[! reduces it to v when exact 
| |D-E| | 


cancellation occurs. This effect is also felt through the 


term ||165-65|!, 


| |D-E | | 

Observe that each term in the braces of (49) has a 
multiplier which is extremely small. That is 6E, 6D, 6S, 
ôq, and اوہ‎ 


roundoff error. The factor ||q|| is generally on the order 


are all of the order of the machine unit 


of unity. From equations (23), (24), (25), and (26) of 


Chapter II, the infinity norm of S is 


20^ 
4 


|١ 


- 0°. (51) 


From equations (l6c), (16d), (21), and (22) of Chapter II, 


the infinity-norm of (D-E) is 
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i + d. 
max le | i 


J 


| |D-E| | 
4 


لا 
De ٦‏ 
í A‏ 


2lu; 
z ue : (52) 


As a consequence of the above facts, the quantity 


T 
] ve! | 6D-déE 
0 D-E 


must be of the order of 1014 for the relative uncertainty to 
be significant. Assuming that ||ó6D-6E||l = 0 and p š v, this 


means that 


v? = 101 4 


OF 


Vs 10? 


before significant roundoff error is propagated to the 
iterate solution. These extremely large condition numbers 
simply do not exist for physical problems of interest in the 


application of the Sy method. 


Estimating the Condition Numbers 

The condition number of the iteration matrix is not 
readily estimated. Although the infinity-norm of (D-E) `!S 
can be calculated precisely, as discussed in Chapter II, 


calculation of || [ወ-5)”፤5]”!| | is not tractable. 
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Consequently an estimate of p is not available nor is an 
upper bound on it. The following discussion assumes that p 
is small enough that its influence on the roundoff error 
bound is less than that of v. 

Estimation of v is possible. Equations (27) and (28) 


in Chapter II reveal that the infinity norm of (5-5) ”? is 


2 R-1 
-] | Tit 5S ٠ (53a) 
| | (D-E) || 2 max4q-! * - + - (Ns 
j j J J J 
Or 
R 
a ہہ"‎ E 
|>-)'|| - Er I 


where R is the number of equi-spaced spatial intervals, A, 





and 
0 < B = | ١ > 1 
፦ i : : 54 
d; (54) 
Here 
1 
_ | AG 
4 ZA (55a) 
and 
t 
2۱۰ + Ag 
_ man 
d: = للا‎ . (55b) 
From (53), 


-ጅ] أ‎ 3 
| | (D-E) |  :ቭፐ።።ዝ +. (56) 
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Consequently from (52) and (53) 


ws 2lujlmax 1-۴ 
Ad; 1-8 (57) 
where { is the index for the minimum absolute ordinate in 
the set 0٢ 


A tight bound on v can be obtained by neglecting کن‎ 





That is, 
2|u; 
0 سا‎ 
Aot 
Since all problems must have Dios E], 
2 
< 
uu (58) 


is a useful estimator for v. For the problems considered 


in Chapter II, ot = 1.0 cm’! and A = 1/3 cm so 
AO e (59) 


Hence (D-E) is a very well-conditioned matrix for the 
problem considered. Consequently roundoff error propagation 
was insignificant in the investigations conducted. In 
practice, roundoff error propagation has not been experi- 
enced by this investigator which renders the Sy algorithm 
unstable. This must be a consequence of the well-conditioned 
state of the matrices (D-E) 7's and (D-E) as well as to the 


roundoff error cancellation property of the iterative method. 
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Effect of Exponential Transform 
The presentation in Chapter III reveals that the 


transform 
IS = TS (60) 


applied to the Sy equations, produces an iterative algo- 
rithm with the same form as that of the unmodified Sy 
algorithm but with a modified source and a different 
iteration matrix, a As a consequence, the roundoff 
error analysis, previously presented, carries over into the 
transform domain in its entirety. That is, expressions 


(25), (37), and (49) are valid for 


l || 
if (D-E) and (D-E) ÍS are replaced by (D-E) and 5 Is 
respectively and q is replaced by its transformed counter- 
part. All of the above quantities are defined in Chapter 
TIT., 
From equations (18) and (19) of Chapter III it is 


clear that 
|| [29-፻] || ። | | (5-5) || (61) 


since terms involving % cancel. From equations (23) and 


(24) of Chapter III, 


R 
- 1 ENEMY | 
D-E - —nrsla ea مس‎ 
| | ) || 27 )62( 
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where 
ሮ/ 
0 > لل - ب‎ > 1 )63( 
0 
and 


2lu.].. -a0%-aalu.].. 


min min 
E 7 ۸ (64a) 


2{u.{_. +40%+aalu .| 


min j ' min 
d; = ہج کہ ہت ےہ‎ . (64b) 


Observe that in problems for which the transform method is 


effective, i.e. ን 0 for all jf, and 


| 1 
< 
d. d. 65 
j 1 D 
of equation (55b) and 


ኘ ሩ 8 (66) 


of equation (54). Since 一 is the predominant factor in 


(62), 
۰ E UA (67) 


Consequently the condition number of (D-E) is smaller than 
the condition number of (D-E). As previously discussed, the 
condition number of this matrix is the predominant factor 

in the stability of the iterative algorithm. Consequently 
it is expected that the transform domain iterative algorithm 


is "better" conditioned than the unmodified Su algorithm. 





CHAPTER VI 
CONCLUSIONS 


The Sy algorithm in one-dimensional plane geometry 
has been formulated within a mathematical framework which 
allows detailed insight into its convergence properties. 
Knowledge of these properties permits the derivation of an 
improved convergence criterion which, if met, guarantees 
that the fractional iterative error is arbitrarily small. 


Utilization of this criterion, 


ር -1‏ دلا 
Yh‏ 


is especially important in problems for predominantly scat- 
tering media for which || ٦١ is close to unity. 

The matrix formulation of the iterative Sy algorithm 
allows calculation of the infinity-norm of the iteration 
matrix, | | (5-5) '5| |, for homogeneous media problems in 
which scattering is isotropic in the Laboratory Coordinate 


system. Imposing the condition 
|| (5-5)”!5|| ሩ ! (2) 


leads to sufficient conditions for which convergence of the 


Sy algorithm is guaranteed. These conditions impose 


|] 





کر 


constraints on the values used for the numerical quantities 


of the problem, A, ul oó, and of, The constraint 


min? 


a < 0 (3) 
O 


is sufficient to guarantee convergence and allows a limited 
amount of negativity in the matrix elements of the iteration 


matrix. The constraint 


A < 0 (4) 
0 
renders the iteration matrix non-negative and is sufficient 
to guarantee convergence for all physically realizable 
problems. The solution for such a system iS non-negative 
for non-negative sources. 


A spatial transform of the type 


۷ x] = Aara (5) 


was used on a special half-range problem. This transform 
places the acceleration parameter, “, in positions within 
the iteration matrix such that each non-zero matrix element 
is reduced. This diminishes the norm of the iteration 
matrix and accelerates convergence of the Sy algorithm. 
Reductions of almost 25% in the number of iterations to con- 
vergence have been realized for sample problems using the 


exponential transform 


(6) َء = سس 
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where a > 0. A practical maximum value of a for this 
transform was found to be 
Ar AS 
17171 E E 
jo man 

A practical useful value of 0 15 much smaller than o* and is 
chosen as a compromise between the user's estimate of 
tolerable discretization error and the reduction in calcu- 
lational effort. It is possible that Sj, problems in higher 
dimensional geometries can be formulated in the mathemati- 
cal context displayed in this investigation for plane 
geometry. If transforms can be found which leave the 
angular discretized SN equations invariant and reduce the 
norm of the iteration matrix, vast savings of computational 
effort can be realized for these problems. Possible trans- 
form candidates are some form of the integrating factor for 
the angular discretized equations describing the problem 
after neglecting the scattering gain term. 

Nothing done here suggests that a transform can be 
found which is useful for all problems in a particular 
geometry. Experience of this investigator indicates that 
it is more likely that each problem in the geometry has a 
"best" transform which can be prescribed based upon the 
Source, boundary conditions, and the prescriber's estimate 
of the form of the solution. 

The spatial discretization error in the S, algorithm 
1S proportional to the second derivative of the solution 


DE The spatial transform 
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AX (8) 


| = UE 


with a > O was used on a one-dimensional problem to decrease 
the discretization error. Discretization errors are shown 
to decrease because the $-domain solution has less curvature 
magni tude than the unmodified Sy p-domain solution. 

The transform used to decrease discretization errors 
must be tailored to each individual problem. The choice of 
a successful transform for a specific problem must be based 
on intuition developed for that problem by experience, 
knowledge of analytic solutions to similar problems, or 
guess. A possible candidate is some variation of the 
integrating factor for the problem. 

An upper bound was derived for the relative uncer- 
tainty in the iterate solution of the Sy algorithm due to 
roundoff errors. This bound is expressed in terms of the 
condition numbers of the matrices (D-E) and ا‎ S. The 
bound reveals that roundoff error Cancellation can occur in 
the iterative process due to the factor |lóp-6E||l. ፲፻ 
cancellation does not occur, however, roundoff errors can 
be magnified during the iterative process if the matrices 
MEE) ۱۹ and (D-E) are extremely ill-conditioned. Such an 
event is improbable. The bound on the fractional uncertainty 
for the iterate solution diminishes as iteration proceeds. 
It does not, however, vanish in the limit of an infinite 


number of iterations. 
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For the problems treated here, (D-E) is well- 
conditioned and its condition number is calculable from the 
problem paramcters. The condition number of the iteration 
matrix is not calculable because || [(5-5]”!5]”!| | is not 
calculable. In the problems treated here, propagation of 


roundoff error was not significant. 
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